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I.  INTRODUCTION 


Earth-orbiting  satellites  are  catalogued  in  the  form  of  orbital  element  sets.  For 
maximum  usefulness  of  these  catalogs  to  militar>’  and  civil  space  systems,  element 
sets  need  to  be  updated  as  frequently  and  as  accurately  as  possible.  A  nonlinear 
filter  was  designed  to  perform  these  updates  in  near-real-time  using  observations 
of  satellites  piercing  a  planar  radar  beam. 

The  filter  predicts  time-varying  orbital  parameters  according  to  a  non-linear 
orbital  dynamics  model.  This  model  includes,  along  with  two-body  Keplerian 
motion,  first  order  oblateness  effects,  both  secular  and  periodic,  simple  decreasing 
exponential  atmospheric  drag  effects,  and  white  noise  on  both  the  dynamic  model 
and  the  measured  observations.  This  simple  model  is  found  sufficient  to  track 
several  selected  satellites,  with  proper  tuning  of  the  filter.  It  should  be  noted  that 
even  this  simple  non-linear  model  gives  rise  to  comparatively  complex  filter 
expressions. 

Applying  this  filter  to  the  orbit  improvement  problem  helps  meet  the  need  for 
continual  determination  of  orbit  decay,  collision  avoidance,  satellite  maneuvers, 
etc.  Current  systems  use  a  daily  batch  least  squares  differential  correction  method. 
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II.  REVIEW  OF  LITERATURE 


This  research  involves  the  integration  of  a  simple,  yet  adequate,  orbital 
dynamics  model  with  a  necessarily  suboptimal  extended  Kalman  estimator.  This 
combination  is  applied  to  a  geometrically  unique  detection  scheme  to  provide  near 
real-time  improvement  of  orbital  elements  for  selected  satellites.  The  accuracy 
level  under  consideration  falls  under  the  umbrella  of  the  general  perturbations 
problem  of  orbital  mechanics,  which  is  covered  by  numerous  textbooks,  a 
sampling  of  which  are  referred  to  here  [1-6].  In  general,  solution  procedures 
involve  a  gravitational  potential  function  (earth  oblateness,  third  body,  etc.) 
describing  its  contribution  to  the  perturbation  of  a  classical  two-body  Keplerian 
dynamic  model.  The  resultant  expressions  take  the  form  of  variations  in  some  set 
of  orbital  parameters.  The  two  broad  divisions  of  orbital  parameter  representation 
in  the  literature  are:  1)  some  set  of  6  orbital  elements  which  describe  the  orbital 
path  as  a  conic  section  and  its  orientation  in  inertial  space,  and  2)  Cartesian 
position  and  velocity  of  the  orbiting  body.  The  6  classical  orbital  elements  were 
chosen  for  their  physical  significance  in  terms  of  how  the  satellite’s  orbit  is  being 
affected  by  various  perturbing  forces.  The  particular  satellite  propagation  theory 
being  used  here  is  loosely  based  on  Brouwer’s  theory  [7-11]  for  a  future 
comparison  basis  with  existing  orbit  improvement  procedures.  The  detection 
scheme  is  an  earth-fixed  planar  radar,  the  crossing  of  which  causes  power  to  be 
reflected  back  to  one  or  more  of  multiple  earth-based  receivers  [12]. 


2 


The  most  poorly  modeled  phenomenon  affecting  low  earth  orbits  is  the 
problem  of  atmospheric  drag.  The  overw'helming  consensus  [13-15]  is  that  even 
the  most  complex  atmospheric  density  models  can  only  be  consistently  accurate  to 
within  about  20%.  The  Kalman  filter  incorporates  a  very  simple  decreasing 
exponential  density  model  [16]  corrupted  by  noise  to  track  drag  penurbations. 

The  extended  Kalman  filter  is  used  in  many  situations  for  nonlinear  parameter 
estimation  [17-19]  and  has  been  applied  to  the  field  of  orbit  estimation  in  the  form 
of  orbital  parameter  improvement  by  differential  correction  techniques  [20,21].  It 
is  known  that  such  a  filter  w'orks  w'ell  W'hen  observational  data  are  plentiful 
[22,23].  However,  proper  application  of  such  a  filter  to  a  highly  nonlinear, 
minimal  observation  scenario,  such  as  that  posed  by  the  fixed  radar  plane  detection 
scheme,  can  also  successfully  "track"  satellites. 
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in.  TWO-BODV  ORBIT  FORMULATION 


Two-body  orbital  motion  for  a  spherical  earth  can  be  described  by  a  set  of  six 
parameters  (three-dimensional,  second  order  problem).  Traditionally,  these 
parameters  take  one  of  two  forms.  One  method  describes  satellite  motion  in  terms 
of  Cartesian  position  and  velocity.  These  vectors  are  propagated  according  to  the 
two-body  non-linear  equations  of  motion.  The  other  takes  advantage  of  the  fact 
that  the  equations  of  motion  describe  an  orbital  path  which  is  a  conic  section, 
which,  for  earth-orbiting  satellites,  is  always  an  ellipse.  Satellite  motion  can  then 
be  analyzed  as  a  mean  angular  displacement  along  that  elliptical  path.  The  second 
method  is  used  here  because  of  the  linear  nature  of  the  mean  angular  motion  and 
its  clarity  of  physical  significance. 

A.  KEPLER  TWO-BODY  DYNAMICS 

The  orbital  path  is  described  as  an  ellipse  of  fixed  shape  and  size  located  in  a 
plane  of  fixed  orientation  with  respect  to  the  earth's  inertial  reference  frame.  The 
quantities  describing  this  path,  the  classical  orbital  elements  [1],  consist  of  five 
constants  and  one  time-dependent  quantity.  Two  of  the  constants  (a,  e)  describe 
size  and  shape  of  the  orbit,  while  the  other  three  (i.  O,  to)  orient  the  orbital  plane  in 
space.  The  time-dependent  quantity  (v)  pinpoints  the  actual  location  of  the 

satellite  on  the  orbital  path.  These  elements  are  (see  fig.  3.1): 
a  (semi-major  axis)  -  defines  size  of  the  elliptical  orbit, 
e  (eccentricity)  -  defines  shape  (non  circularity)  of  orbit. 
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Elliptical  Conic  Section 


.  Classical  orbital  elements. 
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i  (inclination)  -  defines  the  angle  between  equatorial  plane  normal  (K)  and 
orbital  plane  normal  (W) . 

Q  (longitude  of  the  ascending  node)  -  defines  the  angle  between  the  vernal 
equinox  direction  (1)  and  the  point  where  the  satellite  crosses 
the  equatorial  plane  from  southern  to  northern  hemisphere 
(ascending  node). 

CO  (argument  of  perigee)  -  defines  the  angle  between  the  ascending  node  and 
the  satellite's  closest  point  of  approach  (perigee). 

V  (true  anomaly)  -  defines  the  angle  between  perigee  and  the  satellite's  actual 
position. 

The  elliptical  orbital  path,  is  described  by  the  equation  of  a  conic  section 


a(l-e') 

r  = - 

1  +  e  cos  V 


(3.1) 


where  r  is  the  magnitude  of  the  satellite’s  position  vector  at  a  specified  v.  In 
general,  v  is  nonlinear  with  respect  to  time,  so  the  mean  anomaly  M,  w'hich  varies 

linearly  with  respect  to  time  for  purely  Keplerian  motion,  will  be  used  as  a  filter 
state  instead.  M  is  related  to  v  through  a  quantity  called  eccentric  anomaly  E.  It  is 

defined  as  follows  (see  fig.  3.2): 


cos  V  = 


sin  V  = 


E  -  e  sin  E 

(3.2) 

cos  E  -  e 

I  -ecosE 

(3.3) 

sin  eVI-c^ 

1-ecosE 

(3.4) 
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Eqs.  3. 1-3.4  define  the  relationship  between  M  and  r  . 

The  filter  state  vector  will  be 

X,  =[a,  e,  i,  Q,  co,  M.f 

where 


1 - 

M 

i _ 

Ck.i 

Ck 

ik 

f^k.: 

^^k 

f^k.i 

Mk.:. 

_a:%,+M,_ 

(3.5) 


(3.6) 


and  Tk  is  the  time  between  obsen’ations.  For  the  classical  2-body  case,  it  can  be 
seen  from  eq.  3.6  that  state  propagation  is  a  linear  problem. 

In  order  to  improve  the  estimate  of  the  state  vector  using  available 
observations,  it  becomes  necessary  to  compare  the  states  with  the  obser\'ables  in  a 
common  reference  frame.  Satellite  position,  which  is  implicitly  a  function  of  the 
states  through  the  relationship  in  eqs.  3. 1-3.4,  can  be  written 


r 


p 


r  cos  V 
rsin  V 
0 


(3.7) 


where  denotes  a  position  vector  in  the  orbital  (PQW)  reference  frame  (see  fig. 
3.1).  The  remaining  slates  are  used  to  transform  satellite  position  from  orbital  to 
inertial  (UK)  coordinates  using 
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Figure  3.2.  Eccentric  vs.  true  anomaly. 


cojcQ  -  scosQci 

-sox:f2  -  ctosQci 

sQsi 

cwsfi  +  scocQci 

~s(OsQ  +  ccocQci 

-cQsi 

scosi 

ccosi 

ci 

(3.8) 


where  c0  =  cos  0  and  s0  =  sin  0  for  brevity.  (See  App.  A  for  development  of 

C^^) 
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Then  satellite  position  in  inertial  coordinates  will  be 


r  cos  V 
rsin  V 


0 


(3.9) 


Now ,  inertial  satellite  coordinates  need  to  be  related  to  a  measurement  in  the  site- 
based  reference  frame. 


B.  OBSERVATIONS 


The  detection  scheme  modeled  here  is  a  fixed  planar  radar  beam  generated  by 
multiple  continuous  wave  transmitters  along  a  great  circle  path  (see  fig.  3.3).  As  a 
satellite  penetrates  this  beam,  phase  information  is  collected  at  one  or  more  of  the 
receivers.  Antenna  geometn,'  at  each  site  makes  it  possible  to  compute  direction 
cosines,  which  are  the  pseudo-obser\'ations  available  for  element  set  improvement. 
These  obser\'ables  are  (see  fig.  3.4) 


cosa  = 


cosP  = 


tl 

IPl 

PN 

IPl 


(3.10) 


where  East-West  cosine  =  cos  a  and  North-South  cosine  =  cos  P  .  (y  can  be 
calculated  from:  cos^  y-t-cos^  P-i-cos^  a  =  l.) 
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Figure  3.3.  Radar  plane  geometry. 


We  have 


p  cos  a  =  p  •  E 
pcosp  =  pN 


(3.11) 


and 


cos  a 
cosp 


(3.12) 
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Figure  3.4.  Receiver  site  geometry. 


is  the  vector  of  observables  to  be  used  for  orbital  parameter  improvement.  (See 
App.  B  for  development  of  z  in  terms  of  the  orbital  elements.) 

The  position  vector  in  the  site- based  frame  is 

r“=C'>^r'  (3.13) 


where 

cfc(|)  cfstf)  s( 

=  -s(()c(az)-sfc(t)s(az)  c(az)c(})  -  s(az)sf  s(})  s(az)cf  (3.14) 

s(az)s<l)  -  c(az)sfc(()  -s(az)c(j)-c(az)sfs(J)  c(az)cf_ 


and 

£  =  latitude  of  receiver 
(|)  =  Q)e(t-tJ-lon„^ 
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CO-  =  earth’s  rotation  rate 

to  =  time  of  last  Greenwich  crossing  of  vernal  equinox 
lonsite  =  longitude  of  receiver 
(See  fig.  3.5) 

(See  App.  A  for  development  of  C*^). 

Recapitulating, 


z,  =g(x^.T, )  (3.15) 

both  of  which  are,  in  general,  nonlinear  functions  of  the  orbital  elements  \^. 
Assuming  for  now  that  perturbations  to  the  tw'o  body  dynamics  and  observation 
noise  can  be  modeled  as  additive  zero  mean,  white  noise  processes,  the  state  and 
measurement  equations  become 

Xk  =f(x,-nT,.,)  +  w,., 

(3.16) 

Zk  =g(x,,TJ  +  v, 

Kalman  filter  equations  are  then 

g.=p«.,h’[hp«„h'  +  r]' 
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(3.17) 


‘‘a 

p-/;=('-g>h)Pk., 


Figure  3.5.  Receiver  site  geometry. 
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where 


Q  =  E[w,wI] 
R  =  E[v.v:] 


P. 


^^3f(x..T.) 

0X, 

ax. 

(3.]  8) 

'pcosa' 

=  o 

pcospj^ 

(See  App.  B  for  derivations  of  <I>  and  H .) 

In  words,  the  actual  observations,  z^,  are  compared  with  computed 
observations,  ^  (based  on  the  classical  two-body  propagation  model),  resulting 

in  the  Kalman  filter  innovations.  These  are  multiplied  by  the  filter  gain,  which  is 
based  on  the  dynamic  model  and  measurement  uncertainties  and  used  to  produce 
an  improved  estimate  of  the  orbital  elements. 
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C.  RADAR  PLANE  INTERSECTION 

Meaningful  application  of  the  observ'ations  to  orbital  element  set  improvement 
requires  the  ability  to  predict  subsequent  penetrations  of  the  radar  plane  by  the 
satellite  in  question.  A  simple  way  to  do  this  is  to  define  the  satellite’s  position 
vector  in  a  coordinate  frame  referenced  to  the  radar  plane  (See  fig.  3.6).  Then  the 
out  of  plane  component  (Z)  of  position  can  be  checked  for  a  value  of  zero, 
yielding  an  in*plane  condition. 

Satellite  position,  readily  available  in  orbital  frame  coordinates  (r*")  is 
transformed  to  radar  plane  coordinates  by  performing  two  coordinate 
transformations:  one  from  orbital  to  inertial  frame,  the  other  from  inertial  to  radar 
plane  frame,  or 

(3.19) 


where  denotes  a  position  vector  in  the  radar  plane  (XYZ)  reference  frame  and 

se. 


C^  = 


C0_C0 

•p 

-S0 


C0„S0 

>p 

C0 


1' 


0 


-S0„C0  -S0„S0  C0_ 

jp  ip  rp 


(3.20) 


where 

0,p  is  the  inclination  of  the  radar  plane  to  the  equatorial  plane 
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Ion,  is  the  longitude  where  the  maximum  latitude  of  the  radar  plane  occurs, 
which  describes  the  location  of  the  X-axis  of  the  radar  plane 
reference  frame.  (See  fig.  3.6). 

(See  App.  A  for  development  of  C^.) 


Figure  3.6.  Radar  plane  orientation. 


Since  the  only  component  of  satellite  position  of  interest  in  determining  radar 
plane  intersection  is  the  Z-component,  the  expression  can  be  further  simplified  to 
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(3.21) 


r2=[0  0  l]C''^C>Pr'' 

=  [-sG^Pjce  -  se^P4se+ ce^p,  ]r  cos  v 
+[-se  ,p P2C6  -  sO  ^  P5 se  +  cG  ^  Ps  ]r  sin  v 


where 


= 


P  P 


(3.22) 


The  in-plane  condition  of  relevance  to  the  filter  is  that  which  falls  inside  the  field 
of  view  of  one  of  multiple  receiver  sites,  also  placed  along  the  great  circle  path 
circumscribed  by  the  transmitters.  The  lime  epoch  of  this  condition  is  the  time 
close  to  w  hich  we  expect  to  observe  the  satellite  from  one  or  more  of  the  radar 
plane  receivers. 
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IV.  DOMINANT  PERTURBATIONS 


Satellite  orbital  parameters  deviate  from  Keplerian  motion  deterministically 
and  randomly  due  to  many  factors.  Ignoring  for  now  such  intentional  human- 
caused  phenomenon  as  plane  change  and  stationkeeping  maneuvers,  the  largest 
natural  contributor  to  parameter  variations  is  the  nonuniformity  of  the  gravitational 
field  due  to  earth’s  oblateness.  For  low  earth  orbits,  the  next  most  significant 
factor  is  atmospheric  drag. 

A.  GRAVITY  PERTURBATIONS 

The  largest  contributor  to  perturbations  to  the  ^'  .^  lerian  orbital  elements 
arises  from  the  oblateness  of  the  earth  These  perturbations  take  the  form  of  a 
combination  of  secular  and  periodic  variations  of  the  classical  orbital  elements. 
Many  solutions  have  been  developed  for  this  problem,  each  with  its  own  strengths 
and  weaknesses.  All  the  solutions,  however,  are  based  on  the  gravitational 
potential 

U  =  — (4.1) 
r  L  n=2  r  J 

where 

p.®  is  earth's  gravitational  parameter 
r  is  the  distance  of  the  satellite  from  earth's  center 
Jn  are  n*^  order  zonal  harmonics 
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Pn(sin  5)  are  Legendre  polynomials 
6  is  satellite  declination 

and  are  more  or  less  complex  depending  on  the  order  to  which  the  terms  of  the 
expression  are  carried.  This  expression  for  U  already  carries  with  it  a  degree  of 
simplification  in  that  it  assumes  axial  symmetry  for  the  earth.  In  truth,  the  equator 
is  slightly  elliptical,  but  the  main  effect  is  on  the  stability  of  geosynchronous 

satellites,  since  their  orbits  are  in  the  equatorial  plane  [3].  Further  simplification  is 
accomplished  by  dropping  terms  for  n  >  3.  This  is  justified  for  medium  accuracy 

applications  because  of  relative  magnitudes  of  the  Jn  terms 

J,  =1.0826  10'' 

Jj  =-2-10^ 


where  successively  decrease.  So  the  simplified  potential  function  is 


U  =  f3sin^6-l)' 

W  r"  '  2 


(4.3) 


Defining  R  =  13(3  sin ^  6-I)  as  the  perturbing  function,  the  variation  of  that 

function  with  respect  to  each  of  the  orbital  elements  is  [  4  ] 
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^  =  -2r 

3a  a 


aR 


3J. 


[cos  (o(l  +  e  cos  v)(l  -  3  sin M  ■  sin  ^  u) 


de  2(l-e^)r 

-  sin  ^  i  •  sin  2u  •  sin  v(2  -r  e  cos  v)] 

3R  -3 J  2  .  .  .  •  2 

—  =  — —  sin  1  •  cos  1  ■  sin  u 
oi  r 

aR  -3J2  .  2  •  • 

—  =  — ^  sin  1  •  sin  2u 

3(0  2r' 

1^  =  0 

aQ 

aR  _  1  dR 

aM  n  dt 


(4.4) 


To  develop  a  first  order  theory,  the  right-hand  members  of  eq.  4.4  are 
expanded  in  a  Taylor  series  about  the  initial  orbital  element  values,  retaining  only 

the  first  term  of  the  expansion,  yielding  expressions  for  the  derivative  of  each 
orbital  element  with  respect  to  true  anomaly  v.  Integrating  these  expressions 

yields  an  analytical  first  order  solution  in  the  form 

x  =  x„-t-5.x(v)-i-5pX(v)  (4.5) 


where 

X  =  instantaneous  orbital  elements 

Xo  =  initial  orbital  elements 

6sX  =  secular  variations  as  a  function  of  v 

5pX  =  periodic  variations  as  a  function  of  v 
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All  six  orbital  elements  experience  periodic  variations.  Ho\\e\  er.  only  Q.  co 


and  M  experience  secular  perturbations.  These  are 
5.Q  =  -^[cosij2(v-vJ 


6. CO  = 


2p^ 


2--sin^  i, 
2 


(v-vj 


(4.6) 


6  M  =  n , 


1  +  ^ 


(l-3sin'  5J 


% 


(t-tj 


Note  that  8,M  is  expressed  in  terms  of  independent  variable  t,  which  is  the  time  at 
which  true  anomaly  is  equal  to  the  value  of  vin  and  6^(0.  The  periodic 
variations  in  co  and  M  do  carry  l/e„  terms  w’hich  approach  singularity  for  near 
circular  orbits.  This  can  be  explained  by  the  fact  that  co  and  M  are  undefined  for 

circular  orbits  and  ill-defined  for  near  circular  orbits.  Fonunately,  the  combination 
of  CO  -I-  M  obviates  the  singular  condition  by  effectively  forming  one  new  orbital 

element  out  of  co  and  M. 

Long  term  propagation  of  the  orbital  elements,  assuming  a  drag-free 
environment,  is  accomplished  by  applying  only  the  secular  variations,  since 
periodic  variations  will  have  no  net  effect  on  the  elements  over  time.  However, 
the  i>eriodic  variations  will  be  seen  later  to  come  into  play  when  observ'ations  are 
assimilated  into  the  orbital  element  improvement  process.  Inclusion  of  these 
secular  variations  in  the  dynamic  model  yields 
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(4.7) 


where 


f(x,,T,)  = 


3J. 


-7T[cos(iJ](v,.,  -vj 
^Pk 


(0,  + 


3J, 


2p/  L 


2  —  sin* 
2 


(Vk*i  -  vj 


+  n^ 

3 

(l  -3sin^  5^ ) 

ak 

Tk 


(4.8) 


where 

Pk  semi-latus  rectum  of  the  orbit 

Jj  is  the  second  order  zonal  harmonic  (equatorial  bulge). 

The  gravitational  force  model  which  gives  rise  to  these  secular  variations  is 
relatively  simple.  It  assumes  that  the  earth  is  symmetric  about  its  spin  axis  and 
about  the  equatorial  plane.  Noteworthily,  the  first  three  mean  orbital  elements 
experience  no  secular  variation  due  to  earth's  oblateness. 
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It  may  be  possible  to  maintain  identification  of  some  satellites  using  secular 
variations  alone.  Continuous  identification  of  a  satellite,  however,  may  depend 
upon  the  application  of  the  periodic  variations  to  the  mean  elements  at  obser\'ation 
times.  These  periodic  variations  are  assimilated  as  an  additive  (superposition 
assumed)  adjunct  to  the  existing  dynamic  model  and  therefore  are  not  technically  a 
part  of  the  propagation  dynamics.  The  filter  innovations  resulting  from  a 
comparison  of  observables  with  predicted  states  are  then  applied,  via  the  Kalman 
gains,  to  the  mean  elements,  implying  a  prior  removal  of  the  periodic  variations. 
These  periodic  variations,to  first  order,  and  ignoring  terms  of  the  order  of  eo  and 
smaller,  are 


^'^"2p:(l-e^) 


sinM„  •cos[2((o„  +  v)] 


6e 

^  2p^ 


l-|sin'  i„ 

2  2 


cosv+^sin^  i„  •cos(2(0„  +  v) 


7 

+  — sin*  i^  •cos(20i)„  +3v) 


3J. 


Si  p  =  -  ^  sin(2i  „ )  cos[2(a)„  +  v)] 
OPo 


3J, 

S^p 

=  — 7C0si 

p 

4Pp 

3J, 

S^2, 

= - ^COSl 

P 

4p^ 

3J,  f 

SM„ 

P 

4Pol 

(4.9) 
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B.  DRAG  PERTURBATIONS 


Satellites  traveling  in  the  altitude  regime  below  600  km  experience  non- 
negligible  drag  force  due  to  atmospheric  resistance  [3].  This  force  is  well 
represented  as 

FD=ic„-pV=  (4.10> 

2  m 


where 

m=mass  of  satellite 

Cd=2  is  a  dimensionless  drag  coefficient 
A=  effective  cross  sectional  area  of  satellite 
V=  velocity  of  satellite  with  respect  to  atmosphere 
p=  local  atmospheric  density. 


1.  Atmospheric  Drag  Effects 

To  first  order,  drag  force  doesn’t  affect  i  and  Q  [6].  The  main  effects  of 

this  drag  force  are  to  bring  about  secular  decreases  in  a  and  e,  i.e.,  to  decrease  the 
size  and  eccentricity  of  the  orbit,  co  and  M  also  experience  changes  due  to  drag, 

but  these  are  negligible  compared  to  the  oblateness  effects. 

It  turns  out  that 


da 

dM 


m 


l  +  ecosE^^ 


J-ecosEy 


(i-xr 


(4.11) 


where 
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d(l  -e  cosE) 

X  = - 

1  +  e  cos  E 

d  =  — cosi 
n 


n  = 


(4.12) 


and 


=  _  212:1  pa(i-eM 
dM  -  ^  ^  f 


C.A 


m 


(1  +  ecosE) 
(1  -ecosE) 


(1-T)'  cosE 


CdA  -ecosEy^  ,  .  :  c 

— paed  - - -  (l-T)sin  E 

2m  vl  +  ecosEy 


(4.13) 


All  quantities  in  these  two  expressions  are  well  known  (according  to  the  simplified 
model)  except  for  atmospheric  density  p,  the  value  of  which  is  related  to  space  and 

time  in  a  complex  manner.  Atmospheric  models  delivering  as  good  as  15% 
standard  deviation  are  in  the  form  of  tables  or  very  complicated  empirical 

formulae  or  both  [13].  For  simplicity,  it  is  advantageous  to  use  an  analytical 
function  to  model  p.  It  is  assumed  that  variations  in  p  in  the  upper  atmosphere 

(150  km  to  750  km)  are  a  superposition  of  four  individual  variations  [4], 

The  first  variation  is  diurnal,  or,  daily.  This  fluctuation  has  to  do  with 
change  in  solar  flux  intensity  associated  with  moving  from  darkness  to  daylight, 
and  vice  versa.  Density,  depending  on  latitude  and  altitude,  may  change  by  up  to  a 
factor  of  10. 
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The  second  type  is  related  to  the  11  year  solar  cycle  which  is  motivated 
by  intensity  of  sunspot  activity.  Maximum  density,  occurring  ai  the  peak  of 
sunspot  activity,  is  roughly  ten  times  the  mean  density. 

A  third  variation  is  semiannual,  with  density  falling  to  a  minimum  in  July 
and  a  maximum  in  October,  with  a  less  pronounced  minimum  and  maximum  in 
January  and  April,  respectively.  The  range  of  this  variation  is,  generally,  less  than 
1/20  the  range  of  the  1 1  year  solar  cycle. 

The  fourth  type  of  variation  is  nonperiodic  and  unpredictable  and  is  tied 
to  geomagnetic  storms  rising  from  solar  flare  eruptions.  Atmospheric  density  may 
increase  by  a  factor  of  10  from  the  geomagnetic  quiescent  value. 

Over  the  short  term  (days),  barring  magnetic  storms,  the  diurnal  variation 
dominates.  However,  over  several  years,  the  solar  cycle  has  a  more  profound 
effect  on  atmospheric  density. 

The  modeling  approach  followed  here  was  to  start  with  a  very  simple, 
even  crude,  atmospheric  density  model,  including  the  above  considerations  as 
necessary  to  achieve  desired  accuracy. 

2.  Simplified  Drag  Effects 

First,  simpler  expressions  for  changes  in  a  and  e  due  to  drag  can  be 
developed  by  viewing  these  changes  as  a  result  of  orbital  specific  energy  being 
extracted  by  drag  force  [14],  Specific  energy  of  an  orbit,  which  is  conserved  for 
an  unperturbed  Keplerian  orbit,  can  be  expressed  as 


The  change  in  a  for  an  incremenial  change  in  f  is  then 


da  = 


(4.15) 


If  it  is  assumed  that  most  of  the  drag  occurs  near  perigee  (atmospheric  density 
exponentially  decreasing  with  altitude),  then 


from  eq.  3.1,  and 


r  =  rp  =a(l-e) 


V'r  =  pi 


2-1 


=  p(l+e) 


(4.16) 


(4.17) 


from  eqs.  4.13  and  4.15.  For  a  satellite  moving  through  dv  near  perigee  (see  fig. 
4.1),  the  specific  energy  extracted  will  be 


Figure  4.1.  Satellite  moving  through  perigee. 


(4.18J 


d£  = 


F 

--^rdv 

m 


1 


fC^A' 

<  m  . 


pV^rdv 


Substituting  eq  4.17  into  eq.  4.18  yields 


d£  =  _  1  -lpp.(i  +  e)dv 
2v  m  7 


(4.19) 


Substituting  eq.  4.19  into  eq.  4.15  gives  an  expression  for  decrease  in  semi-major 
axis  due  to  a  decrease  in  orbital  energy 


fC 

Ua  = - —  a^(l  -*-e)pdv 

\  m  J 


(4.20) 


Since  the  new  Ci  bit  w  ill  pass  through  the  same  perigee  point,  all  of  da  shows  up  as 
a  decrease  in  apogee  radius,  and 

drp  =  da(l  -  e)  -  ade  =  0  (4.2 1 ) 


by  differentiating  eq.  4.16.  Solving  for  dc  gives 


de  =  —  (1  -e)  =  - 
a 


a(l-e  )pdv 

V  m  ^ 


(4.22) 


3.  Atmospheric  Drag  Characterization 

A  crude  model  for  variation  of  p  starts  with  the  assumption  that 

temperature  and  chemical  composition  of  a  gas  remain  constant  [14].  Then, 
density  is  directly  proportional  to  pressure  according  to 


28 


p  =  kP 


(4.23) 


(Realistically,  k  is  a  function  of  temperature  and  chemical  composition.)  At  an 
altitude  h,  the  decrease  in  pressure  accompanying  an  increase  in  altitude  is 

dP  =  -pgdh  (4.24) 


J— =  -kg|dh  (4.26) 

p.  P  0 

p  =  p 

So,  in  general,  atmospheric  density  decreases  exponentially  with  altitude. 
Eq.  4.26  can  be  rewritten 

P  =  P„e'^"-  (4.27) 

where  Po  O'eference  air  density)  and  ho  (density  scale  height)  change  for  different 
altitude  regimes.  For  example,  using  standard  fixed  values  of  p  from  the  "U.S. 
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Standard  Atmosphere,  1976,"  a  reasonable  fit  can  be  made  by  applying  the  values 
in  Table  4.1  for  the  ranges  of  altitude  shown. 


Table  4.1.  DRAG  PARAMETERS. 


Altitude 

po  (kg/m3) 

Zo  (km) 

80  km<z<120  km 

18.5 

5.8 

120  km<z<140  km 

1.21x  10*3 

11 

1 40  km<z<200  km 

3.00  X  10-6 

21 
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V.  THE  RECURSIVE  ORBITAL  PARAMETER  ESTIMATOR 

An  extended  Kalman  filter,  designed  for  orbital  parameter  estimation,  is  based 
on  the  earth-fixed  radar  plane  detection  scenario  described  in  Chap.  Ill  (Also  see 
fig.  5.1).  The  filter  propagates  the  orbital  elements  in  time  according  to  a  first 
order  nonlinear  dynamic  model  which  takes  into  account  the  dominant 
penurbations  arising  from  earth  oblateness  and  atmospheric  drag.  Upon  receiving 
observations,  in  the  form  of  pairs  of  direction  cosines,  the  filter  calculates  an 
improvement  to  the  current  parameter  estimate.  The  filter  operates  recursively  as 
obseiA'ations  become  available.  The  filter  is  written  in  Matlab  code  and  can  be 
found  in  App.  C. 

A.  FILTER  INITIALIZATION 

The  estimator  is  designed  to  receive  two  data  files.  One  contains  an  initial  set 
of  orbital  elements  along  with  subsequent  pairs  of  available  direction  cosines  (see 
App.  C).  The  other  holds  information  about  the  receiver  sites.  Values  of  pertinent 
constants  are  also  set,  as  are  the  initial  times  necessary'  for  orbit  propagation  and 
estimation  of  observables. 
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Figure  5.1.  Estimator  flow  chart. 
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1.  Input  Data  Files 

Two  input  files  are  required  for  filter  operation.  One  contains  an  initial 
set  of  orbital  parameters,  along  with  subsequent  sets  of  direction  cosines.  This 
information  is  all  tagged  with  corresponding  times  (see  App.  C).  The  last  five 
orbital  elements  occur  in  classical  orbital  element  form,  but,  in  lieu  of  semi-major 
axis,  orbital  p)eriod  is  given.  The  conversion  is  simply 


(5,1) 


where  Tsatis  the  period  of  the  satellite.  Note  here  the  apparent  disappearance  of 
earth's  gravitational  parameter  p®.  This  convenience  is  brought  about  by  the  use 

of  canonical  units  [1],  where 


=1 


DU  I 
TU; 


1  DU 4,=  6378.135  km  and  is  eanh's  mean  equatorial  radius. 

1  TU^=  13.44686  min  and  is  the  time  it  takes  for  a  satellite  traveling  in  a 
circular  orbit  at  a  radius  of  1  DU  to  travel  1  DU  of  arc  length. 

Canonical  units  are  used  in  all  filler  calculations.  The  other  input  file  contains 

information  on  the  receiver  site  locations,a}titudes,  and  azimuths.  The  azimuthal 

measurement  indicates,  for  each  receiver,  the  angular  offset  of  its  coordinate  frame 

from  true  north. 
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2.  Constants 


The  earth's  apparent  rotation  rate  is  initialized  as 

(0^  =  .05883378  ^5  7) 

The  radar  plane's  inclination  to  the  equatorial  plane  and  the  longitude  of  its  center 
line  are,  respectively, 

e„  =33.58310° 

(5.3) 

Ion,  =-101.31348° 

The  radar  plane  is  offset  from  earth's  center  by  about  20  km,  or 

=.0031  DU^  (5.4) 


3.  Time 

All  times  are  given  in  Universal  Coordinated  Time  (UCT)  and  converted 
to  TU®s  for  filter  use.  Tgmwch.which  is  time  elapsed  since  the  Greenwich  meridian 
last  passed  the  geocentric  inertial  I  axis,  is  computed  from  a  reference  angle 
(measured  from  1  to  Greenwich)  at  0  hrs,  1  Jan  92,  of 

6  p.wch  =98.920250931'  (5.5) 

using  the  co^  given  previously. 
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4.  Filter  Matrices 

The  error  covariance  matrix  P  is  initialized  to  reflect  the  level  of 
confidence  placed  in  the  accuracy  of  the  initial  orbital  element  set.  The 
measurement  error  covariance  matrix  R  is  set  to 


(5.6) 


where  =<^cosP  =-0002'  reflects  the  uncertainty  of  the  available  direction 
cosines. 


B.  STATE  PROPAGATION 

The  filter  states  are  propagated  according  to  the  two-body  Keplerian  orbital 
dynamics  perturbed  by  earth  oblateness  and  atmospheric  drag.  The  nonlinear  state 
equations  are 
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and  are  propagated  at  arbitrary  increments  of  •  This  is  continued  until  a 

crossing  of  the  radar  plane  is  detected  by  a  sign  change  of  the  radar-plane-normal 
component  of  satellite  position. 


C.  RADAR  PLANE  INTERSECTION 

Detection  of  a  radar  plane  crossing  calls  a  function  which  iterates  to  an  in¬ 
plane  condition  for  estimated  satellite  position.  This  is  accomplished  using  a 
New'ton-Rapson  iteration  on  the  expression  for  the  plane-normal  component  of 
satellite  position.  The  iteration  calls  for  a  calculation  of  an  approximate  AM 
w'hich  will  drive  Tz  to  near  zero.  This  is  given  by 


AM  = 


(5.8) 


where,  from  eq.  3.21, 


r^  =  X ,  r  cos  v  -f  X  ^r  sin  v 


and 


ax, 

=  — ^rcosv-i-X, 

aM  ' 

ax, 

+  — ^rsin  v-i-X, 
aM  ' 


a(r  cos  v) 
dM 

a(r  sin  v) 
aM 


(5.9) 


(5.10) 


Taking  the  individual  derivatives  yields 
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ax, 

aM 


aX; 

aM 


=  -P,  sin  6_ 

1  rp 

=  -P,  sin 

2  tp 


a(cos6) 

aM 

a(cos6) 

aM 


-P.  sin 

4  rp 

-P,  sin0_ 

5  ip 


a(sin  0) 
aM 

a(sin  0) 
aM 


(5.11) 


where 


and 


where 


a(cos0) 

aM 


=  -0)^a“^(sin  0) 


a(sin  0) 
aM 


=  (O4>a"^(cos0) 


a(rcosv)_  asinE 

aM  1  -  e  cos  E 


a(r sin  v)  _  aVl-e^  cosE 
aM  1-ecosE 


rcos  v  =  a{cosE-e) 
r sin  V  =  aVl-e^  sin E 

aE  ^  1 

aM  1  -  e  cos  E 


(5.12) 


(5.13) 


(5.14) 


This  iteration  is  found  to  converge  to  within  z,  <  e  in  three  to  four  iterations, 
where  e  =  10"‘. 
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The  function  then  determines  whether  the  satellite  is  within  the  field  of  view 
of  one  or  more  receiver  sites  by  calculating  the  in-plane  rx  and  ry  components  of 
satellite  position  and  checking  the  inequality  (see  fig.  5.2) 

-7=^=  >  cos  20'  (5.15) 


It  should  be  noted  here  that  the  field  of  view'  needed  to  be  opened  up  to  pick  up 
some  satellite  observ'ations  by  receivers  at  the  edges  of  the  radar  plane. 


/ 

\ 

/ 

within  field  \ 

/ 

of  view  '  \ 

'X  (Ion  =101”) 

Figure  5.2.  Sensor  field  of  view  . 


D.  ASSIMILATION  OF  OBSERVATIONS 

Once  a  satellite  is  in  the  field  of  view  of  radar  plane  coverage,  the  input  data 
file  is  checked  for  observations  at  or  near  its  time  of  arrival.  If  a  pair  of  direction 
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cosines  is  available,  then  a  measurement  estimate  is  calculated  according  to  the 
development  in  App.  B,  where 


z  = 


(5.16) 


is  also  calculated  according  to  App.  B.  The  plant  noise  covariance  matrix  Q  is 
calculated  using  the  squares  of  expected  levels  of  uncertainty  in  each  of  the  states. 
For  example,  if  the  two-body  problem  is  the  dynamic  model,  then  the  Q  matrix 
would  be 

'o:  0  0  0  0 

0  0  0  0 
0  0  o'  0  0 

Q=  , 

0  0  0  o^  0 

0  0  0  0  o^ 

0  0  0  0  0 

where  o„are  proportional  to  the  magnitude  of  the  largest  unmodeled  pei  turbation, 
in  the  two-body  case,  earth  oblateness.  So,  for  the  two-body  problem. 


0 

0 

0 

0 

0 

^MJ 


(5.17) 
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35, i 
2p^ 

i-1 

2 


2p' 

Now,  the  error  covariance  matrix  can  be  calculated  as 

Py  =^P,.yJ>^-^Q  (5.19) 

A 

H  is  then  calculated  as  given  in  App.  B,  m.iking  possible  the  calculation  of  filter 
gain 

G,  =P^_H"(hP^_H'+r)"  (5.20) 

which  is  then  applied  to  the  residuals  z-z  to  calculate  improved  orbital 

parameters 

It  should  be  noted  that,  throughout  the  filter,  the  quantities  v  and  E  are  needed 
for  various  calculations.  E  is  calculated  fom  M  using  Newton-Rapson  iteration  on 

M=E-esinE  (5.21) 

and  V  is  calculated  by  solving  for 
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then  using 


V  =  cos 


cos  E  -  e 
1-ecosE. 


(5.23) 


to  resolve  the  ambiguity. 


VI.  RESULTS 


This  extended  Kalman  filter  is  designed  to  estimate  the  orbital  parameters  of  a 
satellite,  given  observables  in  the  form  of  a  pair  of  direction  cosines  obtained  upon 
the  satellite's  passage  of  an  earth-fixed,  stationary  planar  radar  beam.  One  day's 
worth  of  data  for  each  of  three  satellites  was  obtained  for  filter  testing  and 
adjustment.  Results  of  the  application  of  three  successively  more  comprehensive 
filter  models  are  compared  with  each  other. 

The  first  filter  models  only  two-body  motion,  the  second  adds  secular 
perturbations  due  to  earth  oblateness,  and  the  third  also  includes  periodic 
perturbations  due  to  oblateness  as  well  as  atmospheric  drag  effects.  Each  filter,  of 
course,  models  unknown  perturbations  as  zero  mean,  white,  Gaussian  noise,  in  the 
tradition  of  Kalman  filtering.  Unfortunately,  none  of  the  available  satellite  data 
include  an  orbit  in  the  drag  regime.  The  filters  will  be  compared  in  the  areas  of 
radar  plane  time  of  arrival,  a  residual  small  sample  statistic,  and  orbital  parameter 
estimation. 

A.  RADAR  PLANE  TIME  OF  ARRIVAL 

With  orbiting  objects  arriving  at  the  radar  plane  on  the  order  of  every  few 
seconds,  it  is  important  that  prediction  of  arrival  times  be  fairly  accurate.  Current 
methods  yield  accuracies  of  within  a  second.  The  best  that  can  be  hoped  for  under 
the  current  detection  configuration  is  about  0.25  second,  which  is  the  sensor  time 
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uncertainty.  After  some  crude  tuning  of  the  three  filters,  a  comparison  of  arrival 
times  for  the  three  satellites  is  shown  in  Table  6.1,  where  At  =  t^j^  ~  ^  • 

It  can  be  seen  that,  in  all  cases  and  for  all  models,  the  first  detection  time  after 
the  estimator  is  started  up  (which  occurs  on  the  order  of  12  hours  after  filter 
initiation)  is  always  mispredicted  by  the  filter  by  a  time  on  the  order  of  seconds. 
How'ever,  it  is  seen  that  the  secular  and  periodic  models  cut  down  that  time  error 
bv  a  factor  of  from  three  to  four.  Subsequent  detections  yield  differences  in  time 
between  calculated  and  observed  plane  intersections  of  less  than  1  second  for  a 
one-orbit-later  detection  and  between  two  and  three  seconds  if  detection  occurs 
next  on  the  "backside"  of  the  orbit  (about  a  half  day  later).  Although  this  is  true 


Table  6,1.  OBSERVED  VS.  ESTIMATED  TIMES  OF  ARRIVAL. 


2-bod y 

Secular 

Periodic 

At (sec) 

1 

Satellite  116 

13.94 

4.29 

4.30 

2.21 

0.23 

-3.64 

-2.61 

-0.63 

-0.32 

0.32 

0.87 

Satellite  117 

15.18 

4.76 

4.78 

-0.10 

0.17 

-0.18 

0.27 

-2.40 

-0.93 

Satellite  118 

3.49 

3.49 

-2.24 

-0.81 

0.32 

0.59 
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for  the  two-body  and  secular  cases,  the  incorporation  of  periodic  perturbations 
yields  subsequent  time  differences  of  less  than  one  second  in  all  cases. 


B.  RESIDUAL  SMALL  SAMPLE  STATISTIC 

The  filter  residual  is  the  difference  between  actual  and  estimated  observations, 
and,  as  such,  is  the  basis  for  a  good  performance  measure.  The  small  sample 
variance  of  the  residual  is  calculated  for  each  case  as  follows: 


G 


2 


(6.1) 


where 

X=t^  (6.2) 

i-i  n 

is  the  small  sample  mean,  Xi  being  the  i*  residual.  Plots  of  these  variances  for  each 
case  are  show'n  in  figs.  6.1  -  6.3. 

It  is  evident  that  the  residuals  for  the  N-S  direction  cosines  are  always  very' 
small.  This  is  because  this  direction  cosine  is  always  close  to  90°  and  the 
geometry  of  the  problem  always  forces  the  estimate  of  this  direction  cosine  to  be 
very  close  to  90°.  However,  the  E-W  (in-plane)  cosines  are  more  interesting  to 
observe.  In  general,  the  trend  is  toward  smaller  residual  variances  from 
observation  to  observation  and  from  simple  to  more  complex  dynamic  model.  The 
only  case  where  this  does  not  hold  true  is  for  Satellite  #117,  where  it  can  be  seen 
that  the  two-body  model  is  sufficiently  lost  after  approximately  eight  hours  to 
yield  higher  residuals  than  at  the  previous  set  of  observations.  Note  that  a  "stack" 
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of  asterisks  denotes  several  observations  gathered  at  the  same  radar  plane  crossing, 
which  are  then  applied  one  at  a  time  for  parameter  improvement. 

It  can  be  concluded  from  this  limited  residual  data  that,  for  the  relatively  short 
times  included  in  the  data,  the  filter  is  performing  acceptably. 

C.  ORBITAL  PARAMETER  ESTIMATION 

The  ultimate  objective  of  this  filter  is  to  be  able  to  predict  what  a  satellite's 
position  will  be  at  some  future  time.  Hopefully,  the  filter  is  predicting  plane 
intersection  times  by  accurately  estimating  orbital  parameters.  Figs.  6.4  -  6.6 
show  time  propagating  values  of  the  four  parameters  that  are  expected  to  be 
varying  most  with  respect  to  their  two-body  propagation  values.  All  three 
satellites  show  similar  results.  In  all  cases,  semi-major  axis  estimates  differ 
greatly  from  model  to  model.  Eccentricity  estimates  for  the  two-body  and  secular 
models  are  very'  similar,  with  the  periodic  model  showing  the  greatest  deviation. 
All  estimates  are  ver>'  similar  for  the  longitude  of  the  ascending  node,  though, 
between  observations,  the  two-body  model  does  not  have  a  mechanism  to 
accurately  propagate  this  element.  Argument  of  perigee  estimates  are  veiy'  close 
between  the  secular  and  periodic  models,  while  the  two-body  makes  very  little 
correction  to  this  parameter. 
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Figure  6.1.  Variance  of  residuals  (sat.  #  116). 
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Figure  6.2.  Variance  of  residuals  (sat.  #  117). 
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Figure  6.3.  Variance  of  residuals  (sat.  #  118). 
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Figure  6.4.  Orbital  element  estimates  (sat.  #  116). 
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Figure  6.5.  Orbital  element  estimates  (sat.  #117). 
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Figure  6.6.  Orbital  element  estimates  (sat.  #  118). 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


Orbital  parameter  estimation  was  attempted  using  an  extended  Kalman  filter 
for  each  of  three  orbital  dynamic  models.  Observations  were  angles  only  and 
made  roughly  twice  a  day  per  satellite  passing  through  the  radar  plane.  It  was 
found  that  the  two-body  orbital  model,  including  secular  and  periodic  variations 
due  to  earth  oblateness,  was  able  to  predict  subsequent  radar  plane  crossings  to 
within  a  second  of  actual  crossings. 

The  filter,  which  includes  a  crude  model  for  atmospheric  drag  effects,  was  not 
tested  on  any  satellites  in  the  drag  regime.  Future  research  should  include  data 
from  satellites  with  altitudes  less  than  600  km.  Also,  longer  term  data  should  be 
obtained  to  check  filter  performance  over  the  long  haul  and  also  to  see  if  the  filter 
could  track  longer  term  disturbances  or  perturbations  to  the  orbit  under 
consideration. 

Tuning  of  the  filter  could  also  be  accomplished  by  obtaining  a  larger  data  base 
of  satellite  orbital  elements  and  corresponding  observations.  This  could  be  done 
manually,  but  it  would  be  more  exciting  to  see  if  an  artificial  neural  network  could 
be  designed  and  trained  to  accomplish  the  same  task. 
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APPENDIX  A 


COORDINATE  TRANSFORMATIONS 

The  calculations  and  algorithms  presented  in  this  paper  necessitate  describing 
vectors  in  several  different  coordinate  frames.  These  frames  are  all  orthogonal. 
Transformation  of  a  vector  from  one  reference  frame's  coordinates  to  another’s  is 
accomplished  using  direction  cosine  matrices  (DCMs). 


Reference  Frames 

The  reference  frames  of  interest  are  as  follows  (see  fig.  A.l): 

Geocentric  inertial  (IJK)  -  The  I  vector  points  inertially  in  the  vernal  equinox 
direction,  with  I  and  J  lying  in  the  equatorial  plane,  and  K 
piercing  the  nonh  pole  (coincident  with  earth’s  spin  axis). 

Orbit-fixed  (PQW)  -  The  orbital  plane  is  the  fundamental  plane,  with  P 
pointing  to  perigee  (point  of  closest  approach),  Q  is  the  semi- 
latus  rectum  direction,  and  W  is  the  orbit-normal  (also  orbital 
angular  momentum  direction). 

Radar  plane  (XYZ)  -  The  radar  plane  is  the  fundamental  plane,  with  X 
pointing  to  the  maximum  north  latitude  point,  Y  lies  in  the 
equatorial  plane,  and  Z  is  the  radar  plane  normal. 

Site-based  (HEN)  -  H  is  local  vertical,  E  is  east,  and  N  is  north.  The  origin 
coincides  with  the  radar  site  location  of  interest. 


Direction  Cosine  Matrices 

Three  DCMs  are  developed  here,  performing  one  orthogonal  rotation  at  a 
time,  then  combining  the  rotations. 
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Figure  A.l.  Reference  frames. 


Orbit-fixed  to  Inertial 

Satellite  position  is  most  simply  expressed  in  orbit-fixed  coordinates, 
but  for  comparison  with  earth-fixed  quantities,  it  is  useful  to  transform  both  to  an 
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intermediate  reference  frame.  The  inertial  coordinate  frame  (in  this  case, 
geocentric  inertial)  is  the  obvious  choice.  Fig.  A.2  shows  the  individual  rotations 
required  to  transform  a  vector  from  orbit-fixed  to  inertial  coordinates. 


Figure  A.2.  Orbit-fixed  to  inertial  rotations. 
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Combining  the  rotations  yields 
=  [rot3][rot2][rotl] 

’  c(OcQ  -  stosQci  ccosQ  +  soxQci  scosi 
=  -soxQ  -  ctosQci  -s(OsQ  +  ctocQci  ccosi 
sQsi  -sicQ  ci 


(A.l) 


Because  the  transformation  is  orthogonal. 


cyp  =[c^ 


T 


ccocQ  -  scosQci 

-s(0cf2  -  ccosQci 

sQsi 

ccosQ  +  scocQci 

-scosQ  -1-  ccocQci 

-sicQ 

scosi 

ccosi 

ci 

P2  P3 
P.  P5  Pe 
Pi  Ps  P9 


(A.2) 


Inertial  to  Site~based 

Fig.  A. 3  shows  the  individual  rotations  required  to  transform  a 
vector’s  coordinates  from  the  inertial  to  a  site-based  reference  frame.  Rotation  #1 
keeps  track  of  earth  rotation.  Rotation  #3  points  out  the  fact  that  the  E  and  N 
vectors  are  not  truly  east-  and  north-pointing.  E  is  in  fact  tangent  to  the  great 
circle  circumscribed  by  the  radar  plane,  and  N  is  orthogonal  to  E  in  the  local 
horizontal  plane. 
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where  0=  to®(t  — to)+ 

(Os,  is  earth's  rotation  rate 

to  is  time  of  last  Greenwich  crossing  of  I  axis 

Ion  is  longitude  of  radar  site 


K 


rot  2  = 


cos-f  0  sin  > 
0  1  0 
-sinf  0  cosf 


where  7  is  latitude  of  radar  site 
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Figure  A.3.  Inertial  to  site-based  rotations. 


Combining  the  rotations  yields 
=  [rot3][rot2][rotl] 

(A.3) 

C^C(t)  sf 

=  -s())c(az)-s^c(l)s(az)  c<()c(az)-s(az)s^s(})  s(az)c/ 

s(|)s(az)  -  s^c4)c(az)  -c(j)s(az)  -  s^s4)c(az)  c(az)c/_ 

Inertial  to  Radar  Plane 

The  transformation  from  inertial  to  radar  plane  coordinates  is  very 
useful  to  the  determination  of  satellite  crossings  of  the  radar  plane.  Rotation  #1 
accounts  for  earth  rotation.  Fig.  A. 4  shows  the  individual  rotations  necessar\  to 
perform  this  transformation. 

Combining  these  rotations  yields 

C0„C0  C0_S0  S0_ 

rp  rp  ip 

c’^  =  -s0  c0  0  (A.4) 

-S0„C0  -S0_S0  C0_ 

rp  ip  rp  J 
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where  6  =  co®(t  -  to )  +  Ion , 

ci>®  is  eanh’s  rotation  rate 

to  is  time  of  last  Greenwich  crossing  of  I  axis 

Ion,  is  West  longitude  of  maximum  North 
latitude  of  radar  plane 


where  6,^  is  inclination  of  radar  plane  to  equatorial  plane 


Figure  A.4.  Inertial  to  radar  plane  rotations. 
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APPENDIX  B 


LINEARIZED  FILTER  MATRICES  6  AND  H 

The  extended  Kalman  filter  can  be  viewed  as  a  predictor-corrector  type 
estimator.  States  are  predicted  using  non-linear  plant  dynamics.  Then  a  correction 
is  computed  by  applying  the  Kalman  gain  to  the  filter  residuals  (observations 
minus  estimates).  Calculation  of  the  Kalman  gain  requires  some  type  of 
linearization  of  the  plant  and  measurement  dynamics  (see  eq.  3.18).  Normally  a 
Taylor  series  expansion  is  used,  keeping  only  first  order  terms  which  are  evaluated 
at  the  best  current  state  estimate.  The  general  form  of  the  linearization  of  a 
nonlinear  discrete  function 


is 


(B.l) 


dx, 

u 

it 

dx. 


it 

0X, 


(B.2) 
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Development  of  Linear  Plant  Matrix  O 

As  will  be  seen,  the  state  matrix  <1>  can  be  as  complex  as  one  lets  it 

become.  Following  is  a  presentation  of  three  progressively  more  complex  models; 

1)  Keplerian  two-body  motion, 

2)  inclusion  of  dominant  first  order  earth  oblateness  gravity  effects,  and 

3)  dominant  atmospheric  drag  effects. 

Two-Body  Problem 

State  dynamics  for  a  Keplerian  2-body  orbit  are  given  as 

(B.3) 

Q, 

CO, 

M,  +a'^*T,_ 

Linearizing  gives 

1  0  0  0  0  0 

0  1  0  0  0  0 

0  0  1  0  0  0 

<i)  = 

0  0  0  1  0  0 

0  0  0  0  1  0 

T  0  0  0  0  1 

L  2 
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Secular  Variations  Due  to  Earth 's  Oblateness 

To  first  order,  the  last  three  orbital  elements  are  affected  by  secular 


variations,  so  that 


=  f(x..Tj 


Q,  -^(cosijAv 

3, r  J  a 

M,+a/2  l  +  -2^(l-3sin'6,)  T, 

L 


(B.5) 


where 


- 

l  +  e^  cos\\ 


sin6^  =sini^  sin{v^ 


(B.6) 
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Linearizing  gives 


1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

6  = 

9Q,., 

9Q,,i 

1 

0 

ao.., 

9a  i, 

9e, 

ai. 

9M, 

9(0, ,, 

9(0 

9(0 

0 

1 

5(0.., 

9a, 

de^ 

9M, 

SM,., 

SM.., 

0 

aM,.: 

aM... 

L  5a. 

9e, 

9i, 

9(0. 

9aM, 

where  the  terms  in  the  fourth  row 

are  found  to  be 

-  =  £l 

cosi. 

ap. 

Av-^ 

aa, 

pI 

k 

9a, 

_  aJj 

cosi. 

AV^ 

9e, 

Pk 

9e, 

9Q,, 

,  3J, 

sin  i , 

Av 

9i, 

2pl 

k 

9Q„ 

.  3J, 

-cosi) 

9v, 

9M, 

2Pk^ 

'  9M, 

(B.7) 


(B.8) 
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where 


3e^ 


-2ae 


(B.9) 


l(l  +  e^)sinE  +  2esinEcosE 
-sin  v^(l-ecosE)^ 

The  fifth  row-  terms  are  found  to  be 

..^Pk 


0COk., 

-  A 

8a  ^ 

■  p: 

8(0,,, 

=  A 

8e, 

Pk 

8a)k*, 

-15J, 

5ik 

2p;  * 

8co,„ 

=  A 

8M, 

2pl 

0a 

^Pk 

8e^ 


dv, 


where 


A 


3* 


i 


(B.IO) 


(B.ll) 
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Sixth  row  terms  can  be  evaluated  by 


aMg,, 

_-9J2 

aag 

2r; 

aMg^, 

_-9J, 

acg 

2r; 

aMg,, 

-9J, 

ai, 

3M.., 

_'9J, 

de,. 


A.  A.  T 


dco, 


A.  sin5,T, 


A.  sin5,T, 


3(sin  8^ ) 

ai, 

a(sin5^ ) 

3(0 


(B.12) 


aM. 


3  aAg  d\’ 

1  +  -Ae  Tg-^— ^ 

2  '  a\\  aMg 


where 


dr,  -2agegA,^ -Pg  cosVg 

aeg  A  7^ 


avg 


-^  +  3sin^6g-^-2rgSin6g 

avg  avg 


a(sin  6g ) 


arg  _  ag(l-eg)eg  sin  Vg 
aVg  (l  +  ecosVg)^ 


d(sin5g) 

avg 


=  sinig  cos(Vg  +0)g) 


(B.13) 
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and 


=  l-3sin^  5^ 
As.  =l-e^ 


A 


6. 


f 

V 


J 


A 7^  =  1  +  e^  cos\\ 


(B.14) 


Development  of  z  and  H 

The  relationship  between  the  measurements  (pseudo-observables)  and 
orbital  parameters  (filter  states)  is  contained  in  the  nonlinear  function 

Zk=g(Xk.T,)  (B.15) 


and  is  the  same  regardless  of  the  state  dynamics  model  being  considered.  The 
obserN'ables  of  interest  here  are  the  north-south  and  east-west  direction  cosines,  or 


cos  a 
cosP 


(B.16) 


and  are  related  to  the  range  vector  p  by 


Pe/ 

cos  a 

II 

/P 

cosp 

Pn/ 

L  /pJ 

(B.17) 
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Note  (from  fig.  B.l)  that  the  HEN  coordinate  system  in  which  p  is  received  is  not 
truly  up,  east  and  north.  H  is  local  vertical,  but  E  is  tangent  to  the  radar  plane 
great  circle,  and  N  is  orthogonal  to  E  in  the  local  horizontal  plane.  This  site- 
dependent  information  is  included  in  the  transformation  matrix  (see  App.  A) 
as  the  angle  az  between  true  east  and  E. 


Figure  6.1.  Site-based  coordinate  system. 


So  Pe,  P^,  and  range  magnitude  p  need  to  be  expressed  in  terms  of  the 
orbital  parameters  in  order  to  form  z  =  g(x,  T).  It  can  be  stated  that 

PhI  rpcosy' 

p=  Pe  =  pcosa  (B.l  8) 

Pn  J  LPcosp. 
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In  order  to  describe  the  direction  cosines  in  terms  of  orbital  parameters,  first 
consider  fig.  B.2.  Earth’s  flattening  has  been  exaggerated  in  the  figure  to  clearly 
show  the  effect  on  radar  plane  orientation.  From  the  figure,  it  is  clear  that 

••  =  r^+R..+P  (B.19) 


where 


rofter  “  20km 


(B.20) 


Figure  B.2.  Relationship  between  satellite  position  and  range. 


R(lat)  is  ihe  latitude  dependent  radius  of  a  point  on  the  reference  geoid  (due  to 
flattening)  [7]  and  is  assumed  to  be  in  the  H  direction,  ronseiis  assumed  to  be 
totally  in  the  negative  N  direction,  and  altsue  is  the  receiver's  altitude.  Combining 
the  elements  of  eqs.  B.18  and  B.19  yields 

’  +PCOSY  ■ 

r“  =  pcosa  (B.21) 

+pcosp_ 

is  the  transformation  of  from  the  orbit-fixed  coordinate  system  to  site-based 
coordinates,  or 

r”  =c’^C^pr'’ 

'H,  H3TP,  J 

=  H,  H3  H,  P,  1 

H,  H,  hJLp,  ] 

r(R7  cos  v  +  Rg  sin  v) 
r(Ri  cos  v  +  Rj  sin  v) 
r(R^  cos  v-»-Rj  sin  v) 

where 

R.  =H,P.-KH3P,-t-H,P7 
R,=H,P,-KH3P3-^H,Pg 
R,  =H7P.-^HgP,-hH,P7 

R,=H,P,+H,P,+H,P.  (B.23) 

R,  =H,P, +H,P. +H,P, 

R,  =H,P,+H,P,+H,P, 
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Combining  eqs.  and  B.22  yields 


'Ph' 

r(R,  cosv  +  Rg  sin  v)-R„^^ 

Pe 

= 

r(Rj  COSV  +  R2  sin  v) 

.Pn. 

r(R ,  cos  V  +  R  j  sin  v)  +  r^g^  _ 

(B.24j 


From  eq.  B.24,  p  can  be  calculated  as 

P  =  VPh+Pe+Pn  (B-25) 

Now,  eq.  B.15  can  be  used  as  the  estimate  of  the  observation  to  compute  filter 
residuals.  But  the  Kalman  gain  depends  on  the  matrix  H,  which  is  derived  by 
linearizing  g(X|^ ,  )  with  respect  to  x^. 

The  form  of  H  is  as  follows 


L  da 


afp> 


dM 

d(W. 


aM 


(B.26) 
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where 


1^_D  — ^ 

p  aa  aa 


da 


1  apn 

1  dp 

p  aM 

p'  dM 

1  apK 

1  dp 

p  da 

P  N  2  ^ 
p  da 

1  ap^.  1  ap 

p  aM  p'  aM 


Following  are  the  details  of  this  development. 

=  (Ri  cosv  +  Rj  sin  vj  — 
aa  da 


ao.,  ,  .  X  ar 

=  (R,  cos  v  +  Rg  sin  v)  — 
aa  da 


^Pn 

aa 


dr 


=  (R.  cosv  +  Rj  sin  v)  — 

da 


da 


ir 


^Ph 

aa 


5Pe 


da  j 


(B.27) 


(B.28) 
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where 


(B.29) 


dr 


l-e‘ 


da  1  +  e  cos  v 


The  expressions  for  the  change  in  p  with  respect  to  eccentricity  e  are  identical  to 
eq.  B.28  except  for  the  replacement  of  with  which  is 


dr  -a(2e  +  cosv  +  e^cosv) 
3e  (1  +  ecosv)^ 


(B.30j 


Since  the  Rn  terms  in  eq  B.24  are  functions  of  the  orbital  elements  i,  Q,  and  to. 
the  expressions  for  the  change  in  p  with  respect  to  these  three  elements  are  as 
follows,  with  the  replacement  of  3^  by  the  appropriate  element: 


^Pe 


dR, 


cos  V  +  ■ 


sin  V 


dPr:  fdR, 


dc 


5; 


cos  V  +  ■ 


dR, 

dg 


sin  V 


^Pn 


dR, 


cos  V  +  ■ 


aRs  .  ^ 

— ^sin  V 

d<;  j 


da  p 


.  ^Ph  ,  _  5Pe  .  ^  ^Pn 
V  dq  dq  d; 


(B.31) 
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The  expressions  for  the  change  in  p  with  respect  to  a  change  in  mean  anomaly 

M  is  more  involved  because  of  the  relationship  between  mean  and  true  anomalies. 
This  relationship  is  repeated  here  for  convenience; 


M  =  E  -  e  sin  E 
cosE-e 


cos  V  = 


1  -  e  cos  E 

sin  Ea/1  -e^ 

sin  V  = - 

1  -  e  cos  E 

Taking  partial  derivatives  as  appropriate  yields 

3pH  3(cosv)  _  3(sinv)1  r  _  , 

— — =  R- — ^  +  Rg — r ^  r +  [R,  cos  v+ R j  sin  v1 


^Ph  _ 

R, 

a(cos  v) 

+  R  g 

a(sin  v) 

aM 

- 

aM 

aM 

app  _ 

a(cos  v) 

+  R, 

a(sin  v) 

aM 

aM 

aM 

^Pn  _ 

a(cos  v) 

+  R, 

a(sin  v) 

aM 

aM 

5 

aM 

r  +  [R,  V  +  R  sin  vl  — ^ 

^  '  ■'aM 


ar 

I  +[R.  cos  v  +  R,  sin  vl - 

J  ^  '  ■'a.M 


(B.32) 

(B.33) 

(B.34) 


(B.351 


where 


9(cosv)  5v  ^E 

- ^  =  -sinv— - — 

dM  an  dM 


a(sinv)  dv  dE 

- =  cos  v  — - — 

dM  aE  aM 


dr 

ae(l  -e^)sin  V  ^v  aE 

aM 

(1  +  ecosv)^  aE  aM 

av  _ 

-(l  +  e  ^ )  sin  E  +  2e  sin  E  cos  E 

aE 

-sin  v(l  -ecosE)^ 

aE  _ 

1 

aM 

1  -  e  cos  E 

(B.36) 


All  the  expressions  for  the  partial  derivatives  of  the  Rn  terms  can  be  found  in 
Table  B.l. 
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Table  B.l.  DEVELOPMENT  OF  TERMS  FOR  H. 


ai 

an 

d(0 

ai 

an 

aco 

aR, 

aR, 

aR, 


=  scosnsiH^  -scocnsiH,  +scociHj 

=  -P,H.+P,H, 

=  P,H,+P,H,+P.H, 

=  scosnsiH ,  -  scocnsiH  ^  +  scociH , 

=  -P,H,+RH, 

=  P,H,+P,H.+P,H, 

=  scosnsiH  j  -  scocnsiH  j  +  scociH  3 
=  -P,H. +P.H3 

=  P3H,+P3H3+P.H3 


aR 


ai 


7^  =  ccosnsiH  ^  -  ccocnsiH  5  +  ccociH , 


an 


=  -P3H,+P3H3 


aco 


=  -P,H,-P,H3-P,H, 


ai 


=  ccosnsiH,  -  ccocnsiH,,  +  ccociH. 


an 


=  -P3H,+P3H, 


aco 


=  -P,H,-P,H,-P,H, 


ai 


ccosnsiH  J  -  ccocnsiH  J  +  ccociH. 


an 


=  -P,H,+P,H3 


2- 

aco 


=  -P.H.-P.H3P,H3 
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APPENDIX  C 


FILTER  CODE  (MATLAB) 

Following  is  the  Matlab  code  for  the  extended  Kalman  filter  described  in  this 
paper.  This  particular  filter  models  the  2-body  orbital  dynamics  problem, 
including  secular  and  periodic  perturbations.  It  requires  2  data  files  for  input  (see 
header  of  filter  code); 

datain  -  contains  a  time-tagged  initial  orbital  element  set  and  an  unspecified 
number  of  time-tagged  pairs  of  direction  cosines  and 
corresponding  receiving  sites 

siteinfo  -  contains  pertinent  information  about  the  set  of  receiver  sites 
comprising  the  radar  plane  detection  system. 

The  filter  calls  3  functions: 

enewton  -  performs  a  Newton-Rapson  iteration  to  obtain  eccentric  anomaly, 
given  a  mean  anomaly  and  eccentricity 

intper  -  performs  a  Newton-Rapson  iteration  to  propagate  satellite  position  to 
an  in-radar-plane  condition. 

twopi  -  forces  angular  measurements  to  a  number  between  0  and  2n. 


Extended  Kalman  Estimator 


The  code  for  these  functions  follows  the  filter  code. 

<^^***lti*)^i:|c*^f^^^^L^^^l^^:^t*********************************7^:^fi************* 


This  estimator  requires  2  input  files 
DATAIN  &  SITEINFO. 


datain=[250  116  82  191626.933 

103.591322  .00756997  66.81264 

116.53101 

83  65024.876  6  14.839  .092 

83  83436.223  5  -76.912  .012 


0  0 

189.71773 

0 

0 


26.77689 


to 


run 
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83  83436.259  4 

83  83436.397  1 

83  83436.408  3 

83  83436.425  6 

83  160057.722  6 

83  174304.847  3 

83  174304.904  5 

83  174304.906  6 

83  174304.908  1 


-63.554  .009  0 

58.962  .062  0 

-57.056  .043  0 

-78.987  .034  0 

64.179  .030  0 

-16.315  .020  0 

-58.495  .007  0 

-62.516  -.001  0 

74.290  .012  0]; 


siteinf 0=  [  1 

2  33.44600 

3  33.33239 

4  33.14672 

5  32.28768 

6  32.04323 


32.57840  -116.97016  1.87473e-5 
-106.99824  2.212449e-4  3.139C0 
-93.55039  8.143107e-6  -4.27563 
-91.02096  6.738l34e-7  -5.66895 
-83.53628  1.118243e-5  -9.71924 
-81.92609  3.866064e-6  -10.5800] 


8.56120 


function [xkeep, Tkeep, Tint, delTint, resid, xint, Pkeep] = . . . 

orbper (datain, siteinfo)  ; 
inpels= [datain (1 ,  : )  ' ; datain (2,  : )  '  ]  ; 
observ=datain ( 3 : length (datain ( ; ,  1 )  )  ,  1 ;  5)  / 

K=inpels  ( 1 ) ; 
for  i=l:7 

lat  (i)=siteinfo(i,2) *pi/180; 
lon(i)=siteinfo(i,3) *pi/180; 
alt(i)=siteinfo(i,4)/ 
az (i)=siteinfo(i,5) *pi/180; 
end 

THIS  PROGRAM  PROPAGATES  A  SATELLITE  ORBIT  AS  A  2-BODY  PROBLEM 
WITH  SECULAR  VARIATIONS  IK  OMEGA,  omega,  and  MEAN  AKOMJkLY .  IT 
CALLS  THE  FOLLOWING  FUNCTIONS: 

TWOPI  -  places  angle  between  0  &  2*pi. 

INTSEC  -  iterates  to  an  in-fence-plane  condition  when  a  fence 
plane  crossing  is  detected. 

ENEWTON  -  solves  for  eccentric  anomaly,  given  values  for  mean 
anomaly  and  eccentricity. 


i^**************************************************************** 

*  DEFINE  EARTH  ROTATION  RATE  * 

S^*-************-^*****-)!******************************************** 

earthrot=. 05883378171654;  %  Define  earth  rotation  rate  (rad/TU) . 

J2=1.082645e-3; 

!^****************************r************************************ 

*  DEFINE  VARIOUS  TERMS  * 

s^************************************************»*************** 

R=[4e-8  0;0  4e-8];  %  Measurement  noise  matrix. 
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Pkkml=eye  (6)  *le-8;  %  Initialize  error  covariance  r.atrix. 

Pkeep ( : , 1 : 6) =Pkkml ; 
resid ( : , 1) = [0; 0] ; 

G=zeros ( 6, 2 ) ;  %  Initialize  Kalman  gain, 

f inc=33 . 58310*pi/180;  %  Define  fence  inclination. 

cosfinc=cos (fine) ; 
sinfinc=sin (fine) ; 

lonx=-10l . 31348*pi/180;  %  Define  location  of  fence  x-axis. 

!^******A*******************»********************»******»******K>t* 

*  INITIALIZE  STATES  (ORBITAL  ELEMENTS)  * 
^***********************************************»»***»»:*»******** 

xnow= [ (inpels (7) / (2»pi* 13. 44689317) ) ^ (2/3)  ; 

inpels (8) ; inpels (9:12) .*pi/180]; 
xkeep ( : , 1 ) =xnow;  %  Save  original  states. 

*  CALCULATE  INITIAL  ECCENTRIC  ANOMALY  * 

O 

EA (1 ) =Enewton (xnow (6) , xnow (6)  ,  xnow (2)  )  ; 

i^****************************************»»**x»»*x*«»»***«******* 

*  FIND  INITIAL  r  and  TA  * 

<^*>t************»*»*»*****************»*»**»*»*ji*itx>t»»***»***»**** 

cosTA= (cos (EA (1) ) -xnow (2) ) / (1 -xnow (2) *cos (EA(1))); 
sinTA=sqrt  ( 1  -xnow  (2)''2)*sin(EA(l))/  (1-xnow  (2 )  *cos  (EA  ( 1 )  )  )  ; 
r=xnow (1) * (1-xnow (2) ^2) / (l+xnow(2) *cosTA) ; 

TA ( 1 ) =twopi (as in (sinTA) ) ; 
if  TA(l)<pi 
if  cosTA<0 

TA(1) =pi-TA(l) ;  %  Perform  quadrant  check 

end  %  for  true  anomaly, 

else 

if  cosTA<0 
TA(l)=3*pi-TA(l)  ; 
end 
e.nd 


%********************-******-************llt**x***»****»***»********** 


*  CALCULATE  TRANSFORMATION  MATRIX  (ORBITAL  TO  INERTIAL)  * 


^******************************l»*****************>r**K»*****»*»***» 


cosperi=cos (xnow (5)  )  ; 
sinperi=sin (xnow (5)  )  ; 
cosascnd=cos (xnow (4) )  ; 
sinascnd=sin (xnow (4 )  )  ; 
costilt=cos (xnow (3) ) ; 
sintilt=sin (xnow (3)  )  ; 


%  Calculate  sines  &  cosines  of 
%  angular  orbital  elements. 


Pl=cosperi*cosascnd-sinperi*sinascnd*costilt; 
P2=-sinperi*cosascnd-cosperi*sinascnd* costilt ; 
P4=cosperi’*sinascnd+sinperi*cosascnd*costilt ; 
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P5  =  -sinper  i  *sinascnci+cosperi  *cosascnd*cost  ilt; 

P7=sinperi *sint ilt ; 

P6=cosperi*sintilt ; 

*  INITIALIZE  TIME  * 

oneJan92=98 . 920250931*pi/180; 

TUperday=107. 067 9334 0974 83/ 
secperTU=806 . 81359/ 

JD=inpels (3) ; 
t ime=inpels (4 )  / 
hrs=fix (time*le-4 )  / 
min=fix  (rem (time, le4) *le-2)  / 
sec=rem (time, le2 )  / 

Tfilter= (hrs*3600+min* 60+sec) /secperTU/ 

Tgnwch=rem (earthrot *TUperday* (Tfilter/TUperday+ . . . 

( JD-1) ) +oneJan92, 2*pi) /earthrot; 

Tcld=Tf liter / 

Tkeep (1) =JD+Tfilter*secperTU/86400/  %  Save  original  time 

^*******************»********************X*»»W**»»»»K*lt5<»»»»*»»** 

*  INITIALIZE  TIME  OF  FIRST  OBSERVATION 

time=observ (1,2)  / 

hrs=fix (time*le-4)  / 

min=f ix  (rem (time, le4) *le-2)  / 

sec=rem (time,  le2)  / 

timieobs=  (hrs* 360 0+min* 60  +  sec)  /  secperTU/ 

^Hlt*ifc*'*^i(r**^r******************iir***<t***************'»**iir'**‘»w»i»'*'*w** 

*  CALCULATE  OUT-OF-FENCE-PLANE  COORDINATE  OF  SATELLITE  POSITION  * 

<^H:***Hr***<r**'A**»***»r***»c*xH:iit**Tk*****'***»'*****»**WT**»ik:iif<t*»»*»»»r*»* 

theta=twopi (earthrot *Tgnwch+lonx)  / 
sintheta=sin (theta)  / 
costheta=cos (theta)  / 

zfence= .0031+ (-sinf inc* cost beta *P1 . . . 

-sinf inc*sintheta*P4 . . . 

+  cosf inc*P7 )  *r*cosTA.  .  .  %  Calculate  current  miagnitude 

+ (-sinf inc*costheta*P2 . . . %  of  out-of-plane  coordinate, 
-sinf inc*sintheta*P5 . . . 

+cosfinc*P8) *r*sinTA/ 

i^*********************************************************»****** 

*  PROPAGATE  ORBITAL  ELEMENTS  TO  NEXT  FENCE-PLANE  INTERSECTION  * 

J2=l .082645e-3/ 

P=l/ 
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n=0 ; 

fcr  k=l:K 

Tloop=pi/8*xnow  (1)  ''1 . 5; 

^*******r*********************»»-«t»*»t**i»***********»**»*»********** 

*  INCREMENT  MEAN  ANOMiALY  &  APPLY  SECULAR  VARIATION 

^★**********************»****»*****»**»*****»**»»**»*K*»********* 

mnmotion=xnow (1 ) ^ (-3/2 )  ; 

xnow (6) =twopi (xnow (6) +mnmotion*Tloop) ; 

EA (k+1) =enewton (xnow (6) , xnow (6) , xnow (2) ) ; 
cosTA= (cos (EA (k+1 ) ) -xnow (2) ) / (1-xnow (2 ) *cos (EA (k+1 ) ) ) ; 
sinTA=sqrt  (1-xnow  (2)  ''2)  *sin  (EA  (k+1 )  )  /  ( 1-xnow  (2  )  *ccs  (EA  (k  +  1 )  )  )  ; 
r=xnow  (1)  *  (1-xnow  (2)  ''2)  /  (1+xnow  (2)  *cosTA)  ; 

TA (k+1 ) =twopi (as in (sinTA) ) ; 
if  TA(k+l)<pi 
if  cosTA<0 

TA  (k+1 )  =pi-TA  (k  +  1 )  ;  %  Perforrr.  quadrant  check 

end  %  for  true  anomaly, 

else 

if  cosTA<0 

TA(k  +  l)  =3*pi-TA(k  +  l)  ; 
end 
end 

»  APPLY  SECULAR  VARIATIONS  TO  omega  &  OMEGA 

^★*****»***»****l«*K****Jr»r***;«*»*»*********»»»K*»»>cx»>t»»**»*»*»»** 

delTAsec=twopi (TA (k+1 ) -TA (k) ) ; 
semilat=xnow (1) * (1-xnow (2) ^2) ; 

Ql=-3*  J2/semilat ''2; 

SOMEGA=Qi *cos (xnow (3) ) 12; 

Somega=  (-Q1)  *  (2-5/2*sin  (xnow  (3)  )  ''2)  12; 

delsOMEGA=SOMEGA*delTAsec; 

delsomega=Somega*delTAsec; 

xnow (4 ) =xnow (4 ) +delsOMEGA; 

xnow (5) =xnow (5) +de Is omega; 

*  *  *  -k  it  -k  -k  *  *  -k  *******  -k  ******  -k  *.  -k  *  *  -k  ***************  -K  *  *  -K  *  *  *  *  -K  *  -K  ******  * 

*  CALCULATE  ATMOSPHERIC  DENSITY  * 

s^*  **************************************************************  * 

k2=3e-6; 
k3=21/6378.135; 
atmdens=k2*exp ( (1-r) /k3)  ; 

^^*  **************************************************************  * 

*  APPLY  DRAG  PERTURBATIONS  TO  a  &  e  * 

^**  *************************************************************  * 
kl=-4*6378.135; 

deldrag= [kl*xnow (1) ^2* (1+xnow (2) ) *atmdens*delTAsec; 

kl^xnow  (1 )  *  (1-xnow  (2) ‘'2)  *atmdens*delTAsec; 

0 ;  0 ;  0 ;  0  ]  ; 
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xn  o  w= xn  o  w-r  de  1  dr  ag ; 

»  CALCULATE  TRANSFORMATION  MATRIX  (ORBITAL  TO  INERTIAL)  * 

!^****»*****************»»*******»**«*****»********K»*K»J<*»jlt»»K*** 

cosperi=cos (xnow (5)  )  ; 
sinperi=sin (xnow (5)  )  ; 

cosascnd=cos (xnow (4 ) ) ;  %  Calculate  sines  &  cosines  of 

sinascnd=sin (xnow (4 ) ) ;  %  angular  orbital  elements. 

costilt=cos (xnow (3) ) ; 
sintilt=sin (xnow (3)  )  ; 

Pl=cosperi*cosascnd-sinperi*sinascnd*costilt; 
P2=-sinperi*cosascnd-cosperi*sinascnd*cost ilt ; 
P4=cosperi*sinascnd+sinperi*cosascnd*cost ilt ; 
P5=-sinperi*sinascnd+cosperi*cosascnd* cost ilt ; 

P7=sinperi *sint ilt ; 

P8=cosperi*sintilt ; 

^**>t*******»***********»**********************K*l»Hl«*»X>t********** 

*  UPDATE  TIME  * 

Tgnwch=Tgnwch+Tloop; 

Tfi lter=Tf ilter+Tloop; 

*  CALCULATE  OUT-OF-FENCE-PLANE  COORDINATE  OF  SATELLITE  POSITION  * 

^***:******xX»**X********X***********************»xxxx***x**»x»*** 

theta (k+1 ) =twopi (earthrot *Tgnwch+lonx) ; 
sintheta=sin (theta (k+1) ) ; 
costheta-cos (theta (k+1) ) / 

2fold=2fence; 

2fence= .0031+ (-sinf inc* cos theta *P1 . . . 

-sinf inc’‘sintheta*P4  .  .  . 

+cosf inc*P7) *r*cosTA. . .  %  Calculate  current  magnitude 
+ (-sinf inc*costheta*P2 . . .  %  of  out-of-plane  coordinate, 
-sinf inc*sintheta*P5 . . . 

+cosfinc*P8) *r*sinTA; 

S^*************************************************X******7k******* 

*  CHECK  FOR  FENCE  PLANE  INTERSECTION  * 

^**************i»:*********lk**************************»************ 

*  NEGATIVE  TO  POSITIVE  ??  * 

^*******************************************x******************** 

if  2fence>0 
if  zfold<0 
n=n+l ; 
negtopos=n; 
xnow=xnow ' ; 

[Tloop, Tf ilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA (k+1) , theta (k+1) , EA (k+l ) ,  z fence]  = . . . 
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inttemp (xnow, r, cosTA, . . . 

sinTA,TA(k4l) , theta (k  +  1) , EA (k  +  1 ) , PI ,  P2 ,  P4  ,  ... 

P5,  P7,  P8,  ‘ignwch,  Tfilter,  Tloop)  ; 
xnow=xnow ' / 
end 
end 

^************************************************»**»»*i»********* 

*  POSITIVE  TO  NEGATIVE  ??  * 

^********************************»******************************* 
if  zfence<0 
if  zfold>0 
n=n+l ; 
postoneg=n; 
xnow=xnow ' ; 

[Tloop, Tfilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA (k+1 ) , theta (k+1 ) , EA (k+1 )  ,  z fence]  =  . . . 
irttemp (xnow, r, cosTA, . . . 

SinTA, TA (k  +  1) , theta (k  +  1) , EA (k  +  1 )  ,  PI ,  P2 ,  P4  ,  ... 

P5, P7, P8, Tgnwch, Tfilter, Tloop)  ; 
xnow=xnow ' ; 
end 
end 

TAold=TA(k+l) ; 

^*llr*********************************************»**************** 

*  SATELLITE  IN  NAVSPASUR  WINDOW  ??  * 

%********************r**»******ilt***************«t*»***»**x*»******* 

if  flag==l 

*  CALCULATE  TIME  FROM  LAST  FILTER  UPDATE  (LAST  OBS)  * 
^************:*^************************************************** 


delobs=t imeobs-Tf ilter; 
while  abs (delobs) <1 
if  length (observ (:, 1 )) >=p 
delobs=t imeobs-Tf ilter; 

delT=delobs; 

delTint (p) =delT*secperTU; 

Tint (p) = (timeobs*secperTU) /3600; 

Tupdate=timeobs-Told; 
if  Tupdate<0 

Tupdate=Tupdate+TUperday; 

end 

Told=timeobs; 
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Tfilter=Tfilter+delT; 

Tgnwch=Tgnwch+clelT; 

!^***»*************************»*********»***»*»x*i«»»»»t»i'»»»jt»»»t»» 

*  INCORPORATE  PERIODIC  VARIATIONS  * 

^*******A**********<t**********»*****************»XK**»**»»***»**X 

xmean=xnow; 

daper=3*  J2*xnow  (1)  /  (2*semilat ''2*  (1-xnow  (2)  '"2)  )  *  .  .  . 

sin (xnow (3))^2*cos(2* (xnow (5) +TA (k+1 ) ) ) ; 
deper=3* J2/ (2 *semilat ^2 ) * ( (l-3/2*sin (xnow (3; ) ^2)  *cosTA+ . . . 
1/4  * sin (xnow ( 3) ) ^2* cos (2*xnow (5) +TA (k+l ) )  +  . . . 
7/12*sin (xnow (3) ) ^2* cos (2*xnow (5) +3*TA (k+l ) ) ) ; 
diper=-3*  J2/  ( 8* semi  1  at'' 2)  *sin  (2*xnow  (3)  )  *  .  .  . 

cos (2* (xnow (5) +TA (k+l ) ) ) ; 
dOMper=3* J2/ (4*semilat'2) *cos (xnow  (3) ) * . . . 
sin (2* (xnow (5) +TA (k+l ) ) ) ; 

dnewper=-3*  J2/  (4*sercilat'2)  *  (l-5/2*sin  (xnov;  (3)  )  '2)  .  .  . 

*sin (2* (xnow (5) +TA  (k  +  l ) ) ) ; 
dtemp=3* J2/ (2*semilat'2) * (1/xnow (2) * . . . 

(l-3/2*sin (xnow (3) ) '2) *sinTA- . . . 

l/4*sin  (xnow  (3)  )  '2*sin(2*xnow(5)  +TA(k+l)  )  -r.  .  . 

7/12*sin  (xnow  (3)  )  '2*sin  (2*xnow  (5)+3*TA(k-rl )))  +  .,. 

1/2* ( (1-3/2* sin (xnow (3) ) '2) *sin (2*TA (k+l ) ) - . . . 

(l-5/2*sin (xnow (3) ) '2) *sin (2* (xnow (5) +TA (k+l ) ) ) + . . . 

3/8* sin (xnow (3) ) '2*sin (2*xnow (5) +4*TA (k+l ) ) ) ) ; 
domper=dtemp; 
dMAper=dnewper-domper; 

varper= [daper; deper ; diper / dOMper ; domper ; dKAper ] ; 
xnow=xnow+varper ;  xnow (6) =twopi (xnow (6))/ 

EA (k+l ) =enewton (xnow ( 6) , xnow (6) , xnow (2 ) ) ; 
cosTA= (cos (EA(k+l) ) -xnow (2) ) / (1-xnow (2) *cos (EA(kTl) ) ) / 
sinTA=sqrt  (l-xnow  (2)  '2)  *sin  (EA  (k  +  l )  )  /  (1-xnow  (2)  'cos  (EA  (k-^1)  )  )  ; 
r=xnow (1) * (1-xnow (2) '2) / (1+xnow (2) *cosTA) ; 

TA (k+l ) =twopi (asin (sinTA) ) ; 
if  TA(k+l)<pi 
if  cosTA<C 

TA  (k  +  l )  =pi-TA  (k  +  l )  ;  %  Perform,  quadrant  check 

end  %  for  true  anomaly, 

else 

if  cosTA<0 

TA(k+l)=3*pi-TA(k+l) ; 
end 
end 

%***★************************************************************ 

*  CALCULATE  ANGULAR  DISPLACEMENT  OF  RECEIVER  SITE  FROM  I -AXIS  * 
%************★***********************★**★**★**★**********★******* 
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phi=twopi (earthrot*Tgnwch+lon (observ (p, 3) ) ) ; 

sinphi=sin (phi ) ; 

sinlat=sin (lat (observ (p,  3) ) )  ; 
coslat=cos (lat (observ (p,  3) ) )  ; 
sinlon=sin (Ion (observ (p,  3)  )  )  ; 
coslon=cos (Ion (observ (p,  3) )  )  ; 
cosperi=cos (xnow (5)  )  ; 
sinperi=sin (xnow (5)  )  ; 
cosascnd=cos (xnow (4 ) )  ; 
sinascnd=sin (xnow (4 )  )  ; 
costilt=cos (xnow (3)  )  ; 
sintilt=sin (xnow (3)  )  ; 

Pl=cosperi*cosascnd-sinperi *sinascnd*cost ilt ; 

P2=-sinperi*cosascnd-cosperi*sinascnd*cost ilt / 

P4=cosperi*sinascnd+sinperi*cosascnd* cost ilt ; 

P5=-sinperi*sinascnd+cosperi*cosascnd*cost ilt / 

P7=sinperi*sintilt / 

P8=cosperi* Sint ilt; 

sinaz=sin (az (observ  (p, 3) ) ) ; 
cosaz=cos (az (observ (p, 3))); 

Hl=coslat *cosphi / 

H2=coslat ’‘sinphi; 

H3=sinlat ; 

H4=-sinphi*cosaz-sinlat *cosphi*sinaz; 
H5=cosaz*cosphi-sinaz*sinlat*sinphi ; 
H6=sinaz* coslat; 

H7=sinaz*sinphi-cosaz*sinlat *cosphi; 
H8=-sinaz*cosphi-cosaz*sinlat*sinphi; 
H9=cosaz* coslat; 

R1=H4*P1+H5*P4+H6*P7; 

R2=H4*P2+H5*P5+H6*P8; 

R4=H7’*P1+H8*P4+H9*P7; 

R5=H7*P2+H8*P5+H9*P8; 

R7=H1*P1+H2*P4+H3*P7; 

R8=H1*P2+H2*P5+H3*P8; 
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*  ESTIMATE  MEASUREMENT  * 

Rsite  =  sqrt  (1- (2-1/296.26)  /296 . 26’'sinlat  ^^2 )  +alt  (cbserv  (p,  3)  )  ; 
rof f set= . 0031 ; 

rhoH=r* {R7*cosTA+R8*sinTA) -Rsite; 

rhoE=r* (R1 *cosTA+R2*sinTA) ; 

rhoN=r* (R4*cosTA+R5*sinTA) +roffset; 

rho=sqrt  (rhoH''2  +  rhoE''2  +  rhoN^2)  ; 

zest ( : , p) = [rhoE/rho; rhoN/rho] ; 

size=xnow (1 ) ; shape=xnow (2) ; tilt=xnow (3) / 
ascnd=xnow (4 ) ; peri=xnow (5) ;MA=xnow { 6) ; 
^***»*******»*-*-*************»********»*****»-****»***»**»**»**»*** 

*  CALCULATE  COVARIANCE  MATRIX  P  * 

J2=l . 082645e-3; 
sinEA=sin (EA (k+1 ) ) ; 
cosEA=cos (EA(k  +  l)  )  ; 
sindecl=sintilt *sin  (TA(k  +  l)  -f-peri)  ; 
sernilat=size’'  (l-shape^2)  ; 

*  CALCULATE  APPROXIMJ^TE  delta  (TRUE  ANOMALY)  * 
^***********»;**********************************»»**»************* 

mnmot  i on=xnow  ( 1 )  ( - 3  / 2 )  ; 

delMA=mnmot ion*Tupdate; 
delTA=delMA; 

Alphi  =  -  (l  +  shape''2)  *sinEA+2*shape*sinEA*cosEA; 

A2phi=l -shape *cosEA; 

A3phi=2-5/2*sintilt^2; 

A4phi  =  l-3*sindecl''2; 

A5phi=l-shape^2; 

A6phi=sqrt (l/size+J2/r^3*A4phi) ; 

A7phi=l+shape*cosTA/ 

dOMEGAda=3*  J2/  (size*semilat ''2)  *cost ilt *delTA; 
dOMEGAde=-6*J2*size*shape/semilat''3*costilt*delTA; 
dOMEGAdi=3*  J2/  (2*semilat''2)  *sintilt*delTA; 

dOMEGAdM=-3* J2/  ^'^*semilat^2)  *costilt *Alphi/  (sinTA*A2phi^3)  ; 

doinegada=-3'  J2/  (size*semilat''2)  *A3phi*delTA; 
domegade=6*  02*  size*  shape /sfeinilat^3*A3phi*delTA; 
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ciomegacii=-15* J2/  (2*semilat''2)  *sintilt*costilt*cielTA; 
ciomegaciM=-3* J2/  (2*seinilat'‘2)  *A3phi*Alphi/  {sinTA’A2phi^3)  ; 

drkda}c=A5phi  /A7phi  ; 

dMAda=3/2*Tupdate*A6phi*  (-l/si2e^2+J2*A4phi’‘  (-3/r''4)  *drkdak)  ; 

drkdek= (-2*size*shape* ( l-shape*cosTA) -semilat *cosTA) /A7phi^2; 
dMAde=3/2’‘Tupdate*A6phi* J2*A4phi*  (-3/r^4)  *drkdek; 

dsindkdi=sin (TA (k+1) +peri) *costilt; 
dMAdi  =  3*Tupdate*A6phi*  {-3*J2/r''3)  *sindecl*dsindkdi  ; 

dsindkdo=sint ilt *cos (TA (k+1) +peri) ; 
dMAdoinega=3*Tupdate*A6phi*  (-3*J2/r''3)  *sindecl*dsindkdo; 
drdTA=semilat*shape*sinTA/A7phi^2; 
dsinddTA=dsindkdo; 

dKiessdTA=-3*  J2/r''2*drdTA-  (-9* J2/r''4*sindecl''2*drdTA+  .  .  . 

6* J2/r^3*sindecl*dsinddTA) ; 
dTAdMA=-Alphi / (sinTA*A2phi^3) ; 
dMAdM= 1 + 3 / 2 * Tupdat e * A6phi * dmes sdTA* dTAdMA; 

dada=l+2*kl*size* (1+shape) *atmdens*delTA; 
dade=kl  *size''2»atmdens*delTA; 
dadM=-kl*size''2*  (1+shape) *atmdens*dTAdMA; 

deda=kl*A5phi*atmdens*delTA; 

dede- 1-2 *kl * size* shape* atmdens*delTA/ 

dedM=-kl*size*A5phi*atmdens*dTAdMA; 

phimat= [dada  dade  000  dadM; 
deda  dede  000  dedM; 

0  0  1  0  0  0/ 

dOMEGAda  dOMEGAde  dOMEGAdi  1  0  dOMEGAdM; 
domegada  domegade  domegadi  0  1  domegadM; 
dMAda  dMAde  dMAdi  0  dMAdomega  dMAdM] ; 

^***************************************************************** 
*  CALCULATE  PLANT  NOISE  * 

^**************************************************************»** 

J2=1.0826e-3; 

semilat=xnow(l)  *  (l-xnow(2)  ''2)  ; 

siga2=  (3*J2*xnow  (1)  /  (2*semilat ^^2*  (1-xnow  (2)  ''2)  )  *  .  .  . 

sin  (xnow  (3)  )  ''2*cos(2*  (xnow  (5)  +TA  (k+1 )  )  )  )  ''2; 
sige2=  (3*J2/  (2* semilat ^“2)  *((l-3/2*sin  (xnow  (3)  )  ^2)  *cosTA+.  .  . 
sigi2= (3* J2/ (8*semilat^2) *sin (2*xnow (3) ) * . . . 

cos (2* (xnow(5) +TA(k+l) ) ) ) ^2; 
mnmotion=xnow  (1) (-3/2)  ; 

MApert=mnmotion*3* J2/  (2*semilat^2)  *  (l-3/2*sin  (xnow  (3)  )  ''2)  *  .  .  . 
sqrt  (1-xnow  (2)  ''2)  ; 

sigOMEG2= (3* J2/ (2*semilat ^2) *mnmotion*cos (xnow (3) ) *Tupdate) ^2; 
sigomeg2=  (3*  J2/  (2*semilat''2)  *mnmotion*  .  .  . 
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(2-5/2*sin  (xnow  (3)  )  ^2)  *Tupdate)  ''2; 
sigKA2=  (MApert ''Tupdate)  ^2; 

Q= [ le-2*siga2  00000 
0  le-3*sige2  0000 
0  0  sigi2  000 
000  sigOMEG2  0  0 
0000  sigomeg2  0 
00000  sigMA2] . *le-3; 

%*************************************************»it************* 

*  PROPAGATE  ERROR  COVARIANCE  * 

Pkk=Pkkml ; 

Pkkml=phimat *Pkk*phimat ' +Q; 

^★**************************^ilr***********i*****'****T*t'********'**n*** 

*  CALL  IN  OBSERVATION  DATA  * 

cosalpha=sin (observ (p, 4) *pi/180) ; 
cosbeta=sin (observ (p, 5) *pi/180)  ; 
z ( : / P) = [cos alpha; cosbeta] ; 

*  CALCULATE  LINEARIZED  MEASUREMENT  MATRIX  K  * 

al=l-xnow (2) ^2; 

a2=l+xnow (2 ) *cosTA; 

a3=2*xnow (2) +cosTA* (1+xnow (2) ^2)  ; 

a 4 =2* xnow (2) *sinEA*cosEA; 

a5=l+xnow  (2)  ‘'2; 

a6=-a5*sinEA+a4 / 

a7=l-xnow (2) *cosEA; 

drda=al /a2 ; 

drde=-xnow  (1)  *a3/a2''2; 
dcosTAdT=-sinTA; 

dTAdEA=  (-a5*sinEA+a4 )  /  (-sinTA*a7''2)  ; 
dEAdMA=l/a7; 

dc  o  s  T AdM= dc  o  s  T AdT  *  dT AdE A  *  dE AdMA ; 
dsinTAdT=cosTA; 

dsinTAdM=dsinTAdT*dTAdEA*dEAdMA; 

drdTA=xnow  (1)  *xnow(2)  *al*sinTA/a2''2; 
dr dMA=drdTA* dTAdEA* dEAdMA; 

dRldi=sinperi*sinascnd*sintilt*H4- . . . 

sinperi*cosascnd*sintilt*H5+ . . . 
sinperi’^costilt*H6; 
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dR2di=cosperi*sinascnd*sintilt*H4- , . . 

cosperi*cosascnd*sint ilt .  .  . 
cosperi*costilt *H6; 
dR4di=sinperi*sinascnd*sintilt*H7- . . . 

sinperi*cosascnd*sintilt »H8+ . . . 
sinperi*costilt*H9; 
dR5di=cosperi*sinascnd*sintilt *H7- . , . 
cosperi*cosascnd*sintilt *H8+ . . . 
cosperi * cost ilt *H9; 
dR7di=sinperi*sinascnd*sintilt *H1- . , . 
sinperi*cosascnd*sintilt *K2-^ .  .  . 
sinperi *costilt  *H3; 
dR8di=cosperi*sinascnd*sintilt *H1- , . . 
cosperi*cosascnd*sintilt *H2+ . . . 
cosperi*costilt*H3; 

dRldOM=-P4*H4+Pl*H5; 

dR2dOM=-P5*H4+P2*H5; 

dR4dOM=-P4*H7+Pl*H8; 

dR5dOM=-P5*H7+P2*H8; 

dR7dOM=-P4*Hl+Pl*H2; 

dR8dOM=-P5*Hl+P2*K2; 

dRldorr.=P2»H4+P5*H5+P8*K6; 

dR2dom=-Pl*H4-P4*H5-P7’*H6; 

dR4dom=P2*H7+P5*H6+P8*H9; 

dR5dom=-Pl*H7-P4*H8-P7*H9; 

dR7dom=P2*Hl+P5*H2+P8’*H3; 

dR8doin=-Pl*Hl-P4*H2-P7’*H3; 

dpHda= (R7*cosTA+R8*sinTA) *drda; 
dpEda= (Rl*cosTA+R2*sinTA) *drda; 
dpNda= (R4*cosTA+R5*sinTA) *drda; 
dpda= (rhoH*dpHda+rhoE*dpEda+rhoN*dpNda) /rho; 

dpHde= (R7*cosTA+R8*sinTA) *drde; 
dpEde= (Rl*cosTA+R2*sinTA) *drde; 
dpNde= (R4*cosTA+R5*sinTA) *drde; 
dpde= (rhoH*dpHde+rhoE*dpEde+rhoN*dpNde) /rho; 

dpHdi= (dR7di*cosTA+dR8di*sinTA) *r; 
dpEdi= (dRldi*cosTA+dR2di*sinTA) *r; 
dpNdi= (dR4di*cosTA+dR5di*sinTA) *r; 
dpdi= (rhoH*dpHdi+rhoE*dpEdi+rhoN*dpNdi) /rho; 

dpHdOM= (dR7dOM*cosTA+dR8dOM*sinTA) *r; 
dpEdOM= (dRldOM*cosTA+dR2dOM*sinTA) *t; 
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dpNciOM=  (dR4dOM*cosTA+dR5dOM’‘sinTA)  *r; 

dpdOM= (rhoH*dpHdOM+rhoE*dpEdOM+rhoN*dpNdOM) /rho; 

dpHdom= {dR7dom*cosTA+dR8dom*sinTA) *r; 
dpEdoin=  (dRldom*cosTA+dR2dom*sinTA)  *r/ 
dpNdom= (dR4dom*cosTA+dR5dom*sinTA) *r; 
dpdom=  (rhoH’*dpKdom+rhoE*dpEdoin+rhoN*dpNdom)  /rho; 

dpHdMA= (R7*dcosTAdM+R8*dsinTAdM) *r+ (R7*cosTA+R8 * 
sinTA) *drdMA; 

dpEdMA=  (Rl*dcosTAdM+R2*dsinTAdM)  *r+  (R1  *cosTA-rR2 * 
sinTA) *drdMA; 

dpNdMA= (R4*dcosTAdM+R5*dsinTAdM) *r+ (R4 *cosTA+R5* 
sinTA) *drdMA; 

dpdMA= (rhoH*dpHdMA+rhoE*dpEdMA+rhoN*dpNdKA) /rho; 

drhoEda=dpEda/rho-rhoE*dpda/rho^2; 

drhoNda=dpNda/rho-rhoN*dpda/rho^2; 

drhoEde=dpEde/rho-rhoE*dpde/rho''2; 

drhoNde=dpNde/rho-rhoN*dpde/rho''2; 

drhoEdi=dpEdi  /rho-rhoE*dpdi/rho''2 ; 
drhoNdi=dpNdi/rho-rhoN*dpdi/rho''2; 

drhoEdOM=dpEdOM/rho-rhoE*dpdOM/rho^2 ; 
drhoNdOM=dpNdOM/rho-rhoN*dpdOM/rho^2; 

drhoEdom=dpEdom/rho-rhoE*dpdoin/rho^2; 
drhoNdom=dpNdom/rho-rhoN*dpdom/rho''2 ; 

drhoEdMA=dpEdMA/rho-rhoE*dpdMA/rho''2 ; 
drhoNdKA=dpNdMA/rho-rhoN*dpdMA/rho^2; 

H=[drhoEda  drhoEde  drhoEdi  drhoEdOM  drhoEdom  drhoEdKA; 
drhoNda  drhoNde  drhoNdi  drhoNdOM  drhoNdom  drhoNdMA] ; 

^★★ilr*************************************************-****-******** 

*  CALCULATE  KALMAN  GAIN  * 

^**************************************************************** 

G=Pkkml *H ' *inv (H*Pkkml *H ' +R) ; 
resid ( : ,p) =z ( : ,p) -zest ( : ,p) ; 


*  CALCULATE  IMPROVED  COVARIANCE  * 
Pkk= (eye (6) -G*H) *Pkkml; 
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H  ir  *  »  » 


Pkeep  (:,p*6  +  l:6*  (ptl)  )  =Pkk.; 

*  CALCULATE  IMPROVED  MEAN  ORBITAL  ELEMENTS  ’ 

S^-K*-k**-k***************1r*****1i*********-K*-k**iif€**-K-K*-^-x-kK*i,*-n-K****** 

xnow=xmean; 

G*resid ( : , p) ; 
xnow=xnow+G*resid ( : , p) ; 
xnow (6) =twopi (xnow (6) )  ; 

xint  (1 : 2, p) =xnow (1:2); 
xint (3 : 6, p) =xnow (3:6) *  180/pi; 

*  CALCULATE  IMPROVED  ECCENTRIC  ANOMALY  * 

^*:***************:******************************»*K«»»»r*»»«>!*»***» 

EA (k+1) =enewton (xnow (6) , xnow (6) , xnow (2) ) ; 

<^**********************************»*************j«*»*******»***** 

*  SOLVE  FOP  IMPROVED  COS  (TA)  ,  sin  (TA)  ,  r,  sin(decl),  r.  * 

S^********************************************!-******************* 

cosTA= (cos (EA(k+l) ) -xnow (2) ) / (1-xnow (2) *cos (EA(k+l) ) ) ; 
sinTA=sqrt  (1-xnow  (2)  ^2)  *sin  (EA  (k+1 )  )  /  (l-xr.cw  (2  )  .  .  . 

*cos (EA(k+l) ) ) ; 

r=xnow  ( 1 )  "  { 1-xnow  (2 )  ''2)  /  ( 1+xnow  (2 )  * ccsTA)  ; 

TA (k+1 ) =twopi (as in (sinTA) ) ; 
if  TA(k+l)<pi 
cosTA<0 

TA  (k+1 )  =pi-TA  (k+1)  ;  %  Perforin  quadrant  check 

end  %  for  true  anorr.aly. 

else 

if  cosTA<0 

TA (k+1) =3*pi-TA (k+1) ; 
end 
end 

TAold=TA (k+1) ; 


i^**********************************************-*-'**************** 


if  length (observ (:, 1 )) >p 
P=P+1; 

tiine=observ  (p,  2)  ; 
hrs=fix (time*le-4 ) ; 
min=fix (rem (time, le4) *le-2) ; 
sec=rem (time, le2) ; 

timeobs= (hrs*3600+min*60+sec) /secperTU; 
del obs=timeobs-T filter; 
else 

delobs=2; 
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end 

end 

end 

f  lag=0 ; 

^**i,**it**-K*-*****-^******1t************-***************-K*yi-K*-K***-t**K* 

*  RESET  TIMES  FOR  NEW  DAY  * 

if  Tgnwch>106 . 7955369153698 
Tgnwch=Tgnwch-l 0  6.7  955369153698; 
end 

if  Tfilter>107 . 0879334 
Tfilter=Tfilter-107 .0879334; 

JD=JD+1; 

end 

xkeep ( : , k+l ) =xnow; 

Tkeep (k+l ) =JD+Tfilter*secperTU/86400; 
end 


In-plane  Iteration 

function [T, Tf ilter, Tgnwch, flag, xnow, r, cosTA, . . . 
sinTA, TA, theta, Enow, zfencej  =  . . . 
inttemp (xnow, r, cosTA, sinTA, TA, theta, Enow, . . . 

PI,  P2,  P4, P5, P7, P8, Tgnwch, Tfil ter, T) ; 
xnow=xnow ' ; 

J2=l . 082645e-3; 

earthrot=. 05883378171654;%Define  earth  rotation  rate (rad/TU) . 

lonx=-101 . 31348*pi/180; 

finc=33.58310*pi/180; 

sinfinc=sin  (fine) ; 

cosfinc=cos (fine) ; 

delT=0; 

zfence=l;  %  INITIALIZE  OUT-OF-FENCE-PLANE  VALUE  OF  POSITION 

!^******************************************************»***** 

%  *  ITERATE  TO  IN-PLANE  CONDITION  * 

while  abs (zfence) >=le-8 

%************ilr******i»**************************************** 

%  *  ADD  IN  PERIODIC  VARIATIONS  * 

i^**^**********************^*************'»**********^********* 

semilat=xnow  (1)  *  (1-xnow  (2)  ''Z) ; 

daper=3* J2*xnow (1) / (2*semilat^2* (1-xnow (2)^2))*... 
sin  (xnow  (3)  )  "'2* cos  (2’^  (xnow  (5)  +TA)  )  ; 
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deper=3*  J2/  (2*semilat ''2 )  *  {  (l-3/2*sin(xnow(3)  )  ''2)  *cc:sTA+  .  .  . 
l/4*sin  (xnow  (3)  )  “'2* cos  (2*xnow  (5)  +TA)  +  .  .  . 

7/12*sin (xnow (3) ) ^2* cos (2*xnow (5) +3*TA) ) ; 
diper=-3* J2/ (8*semilat^2) *sin (2*xnow (3)  )  . . . 

*cos  (2*  (xnow (5) +TA) )  ; 

dOMper=3* J2/ (4  *semilat  ^2 ) *cos (xnow (3) ) *sin (2* (xnow (5 )  +TA)  )  ; 
dnewper=-3* J2/ (4*semilat^2) * (l-5/2*sin (xnow (3)  )  "2)  . . . 

’‘sin  (2’*  (xnow  (5)  +TA)  )  ; 
dtemp=3’*  J2/  (2 ’‘semi  1  at  ^2)  ’‘  (1/xnow  (2)  *  ..  . 

( (l-3/2*sin (xnow (3) ) ^2) *sinTA- . . . 

1/4 ’‘sin  (xnow  (3)  )''2’‘sin(2’‘xnow(5)+TA)  +  .  .  . 
7/12’‘sin(xnow(3)  )  ^2’‘sin  (2 ’‘xnow  (5)+3’‘TA ))  +  ... 
l/2’‘  (  (1-3/2 ’‘sin  (xnow  (3)  )  ^2)  ’‘sin  (2’‘TA)  -  .  .  . 
(l-5/2’‘sin  (xnow  (3)  )  ^2)  ’‘sin(2*  (xnow  (5)  +TA)  )  +  .  .  . 
3/8*sin  (xnow  (3)  )  ^2 ’‘sin  (2*xnow  (5)  +4*TA)  )  )  ; 
domper=dte!np; 
dMAper=dnewper-domper ; 

varper= [daper; deper; diper ; dOMper; domper ; dMAper ]  ; 
xper=xnow+varper ; 

EAper=enewton (xper (6) , xper (6) ,xper (2) ) ; 

cosTAper=  (cos  (EAper)  -xper  (2)  )  /  (1-xper  (2)  ’‘cos  (EAper)  )  ; 

sinTAper=sqrt  (1-xper  (2)  ^2)  ’‘sin  (EAper)  /  (1-xper  (2)  .  ,  . 

“cos (EAper) ) ; 

rper=xper  (1 )  ’‘  (1-xper  (2)  ■'2)  /  (1+xper  (2)  ’‘cosTAper)  ; 

TAper=twopi (asin (sinTAper) ) ; 
if  TAper<pi 
if  cosTAper<0 

TAper=pi-TAper ;  %  Perform^  quadrant  check 

end  %  for  true  anomaly, 

else 

if  cosTAper<0 

TAper=3’‘pi-TAper; 

end 

end 

<^*******************************************»**************** 

%  *  RECALCULATE  PERIODIC  TRANSFORMATION  MATRIX  ’‘ 

$^*************■***************************************■>1****** 

Plper=cos  (xper  (5)  )  ’‘cos  (xper  (4)  )  - .  .  . 

sin  (xper  (5)  )  *sin  (xper  (4) )  ’‘cos  (xper  (3)  )  ; 

P2per=-sin  (xper  (5)  )  ’‘cos  (xper  (4) )  -  .  .  . 

cos  (xper  (5) )  ■‘sin  (xper  (4) )  *cos  (xper  (3)  )  ; 

P4per=cos (xper  (5) ) *sin (xper (4) )  + . . . 

sin (xper (5) ) *cos (xper (4) ) *cos (xper (3) ) ; 

P5per=-sin  (xper  (5)  )  ’‘sin  (xper  (4) )  . .  . 

+COS  (xper  (5)  )  *cos  (xper  (4 )  )  ’‘cos  (xper  (3)  )  ; 


J 
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tAO  H)  rA" 


P7per=sin (xper (5) ) *sin (xper (3) ) ; 

P8per=cos (xper (5) ) *sin (xper (3)); 
stuff=s inf  inc*earthrot* xper  (1 )  ''1 . 5*rper*  .  .  . 
cosTAper* (Plper*sin  (theta)  . . . 

-P4per*cos  (theta)  )  +sinf inc*earthrot *xper  (1)  ''1 .5*  .  .  . 

rper*sinTAper* (P2per*sin (theta) -PSper* . . . 

cos (theta) ) +xper (1) / (1-xper (2) *cos (EAper) ) * . . . 

(-sin (EAper) * (-sinf inc*Plper*cos (theta) - . . . 
sinf inc*P4per*sin (theta) +cosf inc*P7per) + . . . 
sqrt  (1-xper  (2)  ‘'2)  *cos  (EAper)  *  (-sinf inc*P2per*  .  .  . 
cos  (theta) -sinf inc*P5per*sin (theta) +cosf inc^PSper ) ) ; 

zfence= .0031+ (-sinfinc*cos (theta) *Plper- sin fine* . . . 
sin (theta) *P4per+cosf inc*P7per ) *rper* . . . 
cosTAper+ (-sinf inc*cos (theta) *P2per . . . 

-sinfinc*sin (theta) *P5per+cosfinc*P8per) * . . . 
rper*sinTAper ; 
delMA=-zfence/stuf  f  ; 
mnmotion=xnow (1) ^ (-3/2) / 
del2T=delMA/mnmot ion; 

ROP AGATE  MEAN  ORBITAL  ELEMENTS  TO  ITERATIVE  VALUE  OF  TIME  * 

xnow (6) =twopi (xnow (6) +delMA) ; 

Enow=enewton (xnow ( 6) , xnow ( 6) , xnow (2 ) ) ; 

TAold=TA; 

cosTA= (cos (Enow) -xnow (2) ) / (1-xnow (2) *cos (Enow) )  ; 
sinTA=scrt (l-xnow(2) ^2) *sin (Enow) / (l-xnow (2) ... 

*cos (Enow) ) ; 

r=xnow  ( 1 )  *  ( 1  -  xnow  (2 )  ''2)  /  ( 1+xnow  (2 )  *cosTA)  / 

TA=twopi (asin (sinTA) ) ; 
if  TA<pi 
if  cosTA<0 

TA=pi-TA;  %  Perforir,  quadrant  check 

end  %  for  true  anorr.aly. 

else 

if  cosTA<0 
TA=3*pi-TA; 
end 
end 

^★Ht********************************************************** 

%  *  APPLY  SECULAR  VARIATIONS  TO  omega  &  OMEGA  * 

i^************************************************************ 

delTA=TA-TAold; 
if  delTA>pi/2 
delTA=delTA-2*pi; 
end 
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if  delTA<-pi/r’ 
delTA=delTA+2’'pi  ; 
end 

semilat=xnow (1 ) * (l-xnow (2)^2): 

Ql=-3* J2/semilat ^2 ; 

SOMEGA=Ql*cos (xnow (3) ) /2; 

Somega=-Ql* (2-5/2*sin (xnow (3) ) ^2) 12; 

delsOMEGA=SOMEGA*delTA; 

delsomega=Somega*delTA; 

xnow (4 ) =xnow (4 ) +delsOMEGA; 

xnow (5) =xnow (5) +delsomega; 

^*iir***********»******»*****»**************»******»**«**»****»:**** 

*  CALCULATE  ATMOSPHERIC  DENSITY  * 

^*******<r**********************************************»»*****»** 

k2=3e-6; 
k3=21/6378.135; 
atmdens=k2*exp ( (1-r) /k3)  ; 

*  APPLY  DRAG  PERTURPATIONS  TO  a  &  e  * 

!^**********************  **llr*****************»***i.»*»**«>r******** 

kl=-4*6378.135; 

deldrag=  [kl*xr  ^x)''2*  (1+xnow  (2)  )  *atmdens  *0611 
kl *xnwW  ( 1 )  *  (l-xnow  (2)^2)  *atrridens*delT 
0;0;0;C]; 
xnow=xnow+deldrag; 

^*******»!*********************************************»*»***» 

*  RECALCULATE  TRANSFORMATION  MATRIX  (ORBITAL  TO  INERTIAL) 

Pl=cos (xnow (5) ) *cos (xnow (4 ) ) - . . . 

sin (xnow (5) ) *sin (xnow (4) ) *cos (xnow (3) ) ; 

P2=-sin (xnow (5) ) *cos (xnow (4) ) - . . . 

cos  (xnow  (5)  )  ’‘sin  (xnow  (4 )  )  ’‘cos  (xnow  (3 )  )  ; 

P4=cos  (xnow  (5)  )  ’‘sin  (xnow  (4 ) )  + .  .  . 

sin  (xnow  (5)  )  ’‘cos  (xnow  (4 )  )  *cos  (xnow  (3)  )  ; 

P5=-sin  (xnow  (5)  )  ’‘sin  (xnow  (4)  )  + .  .  . 

cos  (xnow  (5)  )  ’‘cos  (xnow  (4) )  ’‘cos  (xnow  (3)  )  ; 

P7=sin  (xnow  (5)  )  ’‘sin  (xnow  (3)  ) ; 

P8=cos  (xnow  (5)  )  ’‘sin  (xnow  (3) )  ; 

^•k  -k  *  -k  *****■»!*  -k  -k  -k  -k  -k  *  *  *  -k  -k  *  -k  *  *  -k  *  -k  ***********  -k  *******************  * 

%  ’‘  IMPLEMENT  TIME  CORRECTION  * 

i^************************************************************ 

delT=delT+del2T; 

theta=twopi (earthrot* (Tgnwch+delT) +lonx) ; 
end 

^^**************************************************  ********** 

%  ’‘  UPDATE  ACTUAL  TIME  OF  INTERSECTION  ’‘ 

%★***★****★************************************************** 
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dP  o»P  dP 


Tgnwch=Tgnwch+delT; 

Tf ilter=Tf ilter+delT; 

T=T+delT; 
xnow=xnow ' ; 

*  CALCULATE  FENCE-PLANE  COORDINATES 

x= (cos fine* cos (theta) *Pl+cos fine* sin  (theta) *P4+  .  .  . 
sinfinc*P7) *r*cosTA+ (cos fine *cos (theta) *F2+, . . 
cosfinc*sin (theta) *P5+sinfinc*P8) *r*sinTA; 
y=  (-sin  (theta) *Pl+cos (theta) *P4) *r*cosTA. . . 

+  (-sin  (theta) *P2+cos (theta) *P5) *r*sinTA; 

^★★**»1l[^*»:Tfc*llr***********^Jkr**********w**********TfcTfc************ 

DETERMINE  WHETHER  ANY  RCVR  MAY  BE  ABLE  TO  OBSERVE  SATELLITE  * 

9'**^***’***********'<t**^***^*************'***»'^********it*»wtttk*'** 

o 

if  x/sqrt  (x^2+y''2)>=cos  (40*pi/l80) 
flag=l  ; 
end 

zfence=0; 


Eccentric  Anomaly  Iteration 

function  [EA] =Enewton (EA, MA, e) 

if  MxA<EA 

if  EA>1 . 5*pi 
MA=MiA+2*pi ; 
end 
end 

errE=l ; 

while  abs (errE) >=le-10 
errE=EA-e*sin (EA) -MA; 
delEA= (-1/ (l-e*cos  (EA) ) ) *errE; 
EA=twopi (EA+delEA) / 

MA=twopi (MA) ; 
end 
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Anpj*  Between  0  &  27C 


%  This  function  will  take  any  angle  (in  radians)  and 
calculate 

%  its  equivalent  between  zero  and  2*pi 

function  phi=twopi (phi) 

if  phi>=2*pi 

while  phi>=2*pi 
phi=phi-2*pi; 
end 
else 

while  phiiO 
phi=phi+2*pi/ 
end 
end 
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